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TWO-DIMENSIONAL PLASTIC STRESS SYSTEMS WITH 
ISOMETRIC PRINCIPAL STRESS TRAJECTORIES 
By P. F. NEMENYI and A. VAN TUYL 
(U.S. Naval Research Laboratory, Washington, D.C. 

U.S. Naval Ordnance Laboratory, White Oak, Silver Spring, Md.) 
Received 24 April 1951] 


SUMMARY 

Carathéodory and Schmidt, in their paper on the maximum shear trajectories in 
plane plastic stress with the yield condition 7 —= (¢,—0,)/2 — constant, discussed the 
possibility of these trajectories forming an isometric net. The authors generalize the 
above study to an arbitrary yield condition f(o,,0,) = 0. A complete discussion 
shows that only for three distinct families of yield conditions, depending on from 
two to four parameters, do there exist isometric nets satisfying the conditions of the 
problem beyond the trivial nets consisting of concentric circles and concurrent 
straight lines. Among these yield conditions the von Mises parabolic yield condition 
is found to be included. The linear yield condition is discussed in detail; it is shown 
that + o+constant and o constant are degenerate cases, the latter falling 
under elastostatics. The complex functions Z(z) characterizing the possible isometric 
trajectory patterns are given either explicitly or by equations involving quadratures. 
It is shown that among the isometric nets which are possible for certain yield 
conditions are all those which go into themselves either under a rotation or under 
a translation. 


1. Introduction 


CARATHEODORY and Schmidt (1) investigated the maximum shear trajec- 
tories in plane stress for which + = (o,—o,)/2 is everywhere constant. 
Among other things, they found the complete family of isometric nets which 
can serve as maximum shear trajectories in such a plane system. The prin- 
cipal stress trajectories are, of course, the 45° trajectories of these curves. 

Both in plane strain and in plane stress systems the equation 7 = con- 
stant is consistent with a wide variety of spatial \ield conditions in addition 
to the simplest yield condition maximum sheariig stress equals constant. 
The latter goes back to Coulomb, and has been widely used by Prandtl, 


Nadai. Trefftz, Hencky, and others. A yield condition involving only o, 


and o,, as, for example, 7 = constant, exists in general for a system of 
stresses in a three-dimensional body whenever there is a second relation 
G(o,, 5,63) = 0 connecting the principal stresses in addition to the spatial 
yield condition F(o,,,,0,) = 0; the two-dimensional yield equation is 
then the result of eliminating o, between F = 0 and G = 0. It is easily 
shown that the relation 20,—o,—o, = 0 holds in the case of plane strain 
at right angles to the direction of o, in an incompressible material (2, p. 331, 
for example), while in the case of plane stress, o, = 0 by definition. Thus, 
(Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 1 (1952)] 


5092.17 


B 
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7t = constant, is consistent with any spatial yield condition of the form 
F(o,—o,,0,—o3) = 0 in plane strain, and with any yield condition 
F(o,—o,,03) = 0 in plane stress. 

The von Mises constant strain-energy condition, represented by the 


9 9 


cylinder (of radius k,/%), (¢;—o)?—(o,—99)(0,—o3)+(0,—@3) k?, gives 
37? = k* in the case of plane strain. In plane stress, on the other hand, it 
leads to the von Mises elliptic yield condition, of—o,0,+03 = k*?. On 
introducing the variables o (o,-+-0,)/2 and 7, this becomes 

T + 4/{(k?—o?)/3}. 
It was suggested recently by von Mises that the latter may be replaced 


for convenience by a parabolic yield condition of the form 


9 9 l 
+ | q2g2 : 
‘ 4a? 


The spatial yield condition of St.-Venant, both in plane strain and in plane 


stress, gives a two-dimensional yield condition which is represented by a 
hexagon having one pair of sides parallel to the o-axis. It is therefore, 
consistent with 7 = constant in both of these cases for a limited range of oc. 

The present writers, having given a new demonstration of the result of 
Carathéodory and Schmidt in terms of principal stress trajectories, pro- 
ceeded to find all isometric nets consistent with a general linear yield 
function, and with the two yield conditions of von Mises just mentioned. 
However, they found that a discussion for an arbitrary two-dimensional 
yield equation 7 = f(c) is more profitable, and succeeded in finding all 
yield functions f(c) with which isometric principal stress trajectories 
beyond the trivial ones are consistent. In the course of this analysis, all 
more special problems proposed originally find their solutions. 

If c, and o, depend only on x and y, they evidently satisfy the same 
equations whatever o, may be. The problem which we state in the title 
in terms of plane stress is therefore mathematically the same for all 
relations G(o,,¢, 0,3) = 0, hence identical for plane strain and plane stress. 
However, the most interesting yield equations have simple physical inter- 


pretations only in the case of plane stress. 


2. General equations of the problem 
Referred to the principal stress trajectories as coordinates, Lamé’s 


equations of equilibrium take the form 


l l ¢ lélogr 1 Go 
(r—o) . ——s (L,) 
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TWO-DIMENSIONAL PLASTIC STRESS SYSTEMS 3 
where the directions of s, and s, in that order form a right-handed system, 
and p, and p, are radii of curvature. For the isometric net characterized 
by the complex function Z(z) = $+, associating ¢ with the o,-trajec- 


tories and writing |\dZ/dz V, the well-known relations 
l Clog V | élog V 
Py 08s Pe o8) 


are valid. Substituting in (L,) and (L,), we obtain the simultaneous 
differential equations 


Clog V2 c log tr l do 
O85 O85 T 08,” 
Clog V2 Clogr , 1 da 
O8; 08; T 08, 
Integrating and substituting s = f(c), we find 
: ’ ~< 
log V2 log f(c) | . ca _ 2G), 
J S(e) 
: 1 : 
log V2 log f(c)- a + 2F(¢). 
J S(e) 
Finally, log V* = log f(a) + F(¢)— G() (1) 
* do : ; 
F(d) G(ud). (2) 
(a) 


3. The trivial solutions 

The relations (1) and (2) can be satisfied for any yield equation if either 
F or Gis constant. The two cases are essentially identical, so let us assume 
that G is constant. Then from (2) it follows that o is a function of ¢ alone, 
and by virtue of (1), log V? = a¢d-+-b. Hence, for every yield function f(c), 
any Cartesian net and any net consisting of concentric circles and con- 
current straight lines satisfy the conditions of the problem, and an appro- 
priate o-distribution can be found which is consistent with the equations 
ofequilibrium. For the Cartesian nets, o is obviously constant throughout 
the field. For the net consisting of concentric circles and concurrent lines, 
it follows immediately from (L,) and (L,) that the o-distribution is given 


by the equation 
; P f’(c) t 


I(e) 


2 log? 


l 
do. 


lhe minus sign is valid when we identify the concentric circles with the 
+,-trajectories. 
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4. General discussion 


On denoting the integral on the left side of (2) by u(c), (1) can be written 
in the form 


log V? = logflu-"F+G)|+ F—G = A(F+G)4+ F-G 


when we choose a definite branch of the (at least double-valued) function 
f(c) in a suitable range of co. Both f(c) and o (also F and G) must have 
first derivatives in order to satisfy the equations of equilibrium, hence w(c) 
can be inverted in the neighbourhood of almost every value of o. The 
problem thus consists first in finding all harmonic functions of the form 
A(F+-G)+ F—G, then finding the corresponding yield function for each 
from equations (1) and (2). 

It follows readily that f(c) is analytic when the principal stress trajec- 
tories form an isometric net. Differentiating A(#+G)+ F—G separately 
with respect to ¢d and ys, we find that | A’( f+ G)+-1|F’ and|[A’(F+G)—1]@’ 
are both harmonic, and therefore analytic. The first shows that A’(/'+ @) 





is an analytic function of ys, the second that it is an analytic function of ¢. 
Hence F, G, and A’(F+) must be analytic separately, and so A(x) is 
analytic. From (2) it follows that f(c) is analytic. We do not consider 
whether different isometric patterns can coexist in adjacent regions when 
the yield function is piecewise analytic. 
If A(F+G)+ F—G is harmonic, the differential relation 
A(F+¢@)= 2% Fete aria (D) 
f 2 (’2 fk 2 (2 

must be satisfied. In the investigation which follows, we first find all 
combinations A(u), F, and G for which (D) is satisfied. Having found A (uw), 
we then determine f(c) by the equations 


ota { exp A(u) du. (f,) 

| du , 

(f,) 

J}(o) do x 

which we obtain from A(u) = log f(¢) and u | do/f(c). Finally, if 


A(F+G)+F—G is the real part of 29(Z), the most general complex 
function Z(z) belonging to the yield function f(c) is obtained from the 
equation ¥ 

-a,+id, (c,+-0¢y) | exp| —g(Z+-b,+-ib,)] dZ. 


The constants b, and 6, affect only the parametrization of the curves 
d(x, y) = constant and u(a,y) = constant; a change in the constants 4, 


and a, merely translates the curves, while a change in c, and c, rotates 
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TWO-DIMENSIONAL PLASTIC STRESS SYSTEMS 
and stretches them about the origin. We can therefore omit all six con- 
stants without essential loss of generality. 


5. The non-trivial solutions 

The right side of (D) is seen to be a function of F+G alone. In the 
following discussion of (D) we distinguish two main cases : (I) the two terms 
on the right are functions of /+G separately, and (II) they are not. 
The symbol J(f), for a given function f(¢,%), will always denote the 


Jacobian , 7 
P Leo 
ibe 
db, obs 


In the first case mentioned, the relations 


( ASE ] 


Fr" —G PF’ +.q” 
J 0. J 0 (3) 
F2G? PF? G? 
hold simultaneously. By adding and subtracting, we find 
i yy" ( Pad 
J\ I 0. P 0. 
P22] PF? G? 
That is, G'S FR" (F’2+- G’*)—2F' F"(F"—G’"): 0) (4) 
and FSG" (F’2+-G’2)+.20'@"(F"—G’y: = 0. (5) 
It follows from (4) that (F G”) | (F’24-G’?) is a functien of d alone, and 


from (5), that it is a function of % alone. Thus a necessary condition for 
3) to hold with real functions F and G is F”—G"+a(F’+ G2) 0. a 

| being a real constant. This is equivalent to the pair of equations 
F"+aF"? = 6, (6) 


G aG’2 bh. (7) 


We now complete the discussion of Case I by investigating separately 
the four possibilities a = b = 0;a = 0,b 40;a 40,6 = 0; anda + 0, 


b 0, each of which leads to a non-trivial family of solutions of (D). 


| I 1): a b () 
We have F A,¢+B,, G A,w+ Bb,. Hence A”(u) 0, and 


From (f,) and (f,) we immediately find f(c) vo--constant. The trajec 
| tories in this case are obviously logarithmic spirals. 
([2): a 0. b 0) 
The general solutions of (6) and (7) are F sb(¢-+-a,)*+-constant and 
G $b(%-+-a,)?+-constant. These give A” A’/u, which has the general 


solution 


A(u) vlog w+-B. 
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From (f,) and (f,) we find the yield function f(o) = A,(o+A,)*(*+», 
Omitting unessential constants, the corresponding complex function is 
obtained from the equation 

Z | Z-*exp($bZ?) dZ. (8) 
Setting « = 0 gives tr = f(c) = constant, which is the yield relation used 
by Carathéodory and Schmidt. The above equation with « set equal to 
zero gives a family of nets corresponding to 7 = constant. 
(I 3) and (I 4) 


For a + 0, b = 0, we find 


: | , | 
f : log(?-++-a,)+constant and G ae log(s-+-a,)-+-constant. 
Fora ~0,6 +0, 


, l 
i logeosh ,/(ab)(¢-+-a,) + constant 
a 


. 


| | 
while G —~ logeos y(ab)(ip 1-5) 4 


constant. 


The function A(w) corresponding to the second case does not contain b 
at all; its differential equation is the same as that for the first case with 
u replaced by u+(iz/2a). The form of A(w) is therefore the same in both 
cases, but the arbitrary constants occurring must be restricted differently 
in order for A(F#+G)+ F—G to be real. We have 


I ; , 
A(u) log} (e*"”— 1)!+8(ee" 4 1)1-4\__ y+ constant, (9) 
a 


where & is purely imaginary in the first case and real in the second. The 
corresponding yield function is given by the simultaneous relations 


ota C ( (eau | a B)/a(gau ‘2 )a B)/ag—u du. 
l du 
. = - , (10) 
f(c) do 


For any yield function defined by this pair of equations with B purely 
imaginary, there is a family of isometric nets corresponding to b = 0; 
with B real, there is a family of nets corresponding to b 4 0. These cases 
coincide for B = 0, so two distinct families of isometric nets belong to 
these yield functions. It is seen that the two choices B = 1, a = ay and 
B = 0, a = a,/2 give the same function A(w) (and therefore the same yield 
function f(c)), hence for B 1 also, any yield function is consistent with 
two distinct families of isometric nets. It is easily shown that for no other 
value of B can A(u) coincide with an A(u) for B = 0. 
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Writing B = fi, the function g(Z) corresponding to (13) is equal to 
(1/a+iB/a)log Z+-constant; excluding the case a 1,8 = 0, and including 
unessential constants for the moment, Z(z) is given by the equation 


(b, | ib.) ( Z-A+Bi/a qd Z (c,+ ic,)Z) (1+Bij/a, (11) 
When a 1, f= 0, 


z+a,+ia, logZ. or Z (c, +-tc,)e*. (12) 


Because of the arbitrary complex factors in (11) and (12), the isometric 


nets defined there can serve either as principal stress or as maximum shear 
trajectories. 


In the second case, b ~ 0, B must be a real constant 8. Replacing the 
arbitrary constant ,/(ab)/2 by c, we obtain (omitting unessential constants) 
q(Z) log{sinh"+9)4(¢Z)ecosh"-#)'4(eZ)}, 


and Z(z) is found from the equation 


| sinh-“+4(¢Z)cosh®-Y/4(cZ) dZ. (13) 
It is evident that the choice B lla 1, ¢ L/m?, x 1/2m? 


in (10) gives the parabolic yield function of von Mises, f(a¢) = m*o?— (1/4m?). 
Hence, two distinct families of isometric nets are consistent with this yield 
function. 


Similarly, substituting B= 0, a 1, we obtain the yield function 


f(c) , {(o+-a«)?—4C?}, represented by a hyperbola.t+ However, we can 


show that the elliptic yield function of von Mises, f(o) = ,/{(k?—o?)/3}, 
does not belong to the preceding family, and is therefore not consistent 
with any isometric nets beyond the trivial ones. Indeed, differentiating 
both sides of (f,) with respect to u, we see that the left side is periodic 
and real with a real period for the von Mises yield function, while the 
right side can be seen not to satisfy these conditions for any choice of the 
parameters. 


Case Il 

We assume that neither of the expressions (F” —G")/(F’?+G’?) and 
(F"+-G")/(F?+G") is a function of F+G. If the first is a function of 
F+G but not the second, it follows from (D) that A’(w) = 0, while if the 
second is a function of + G, the first must be a function of F+G neces- 
sarily. 

We show in the following that only a very restricted class of yield 


It seems unlikely that there exists any material having a yield function which is 
monotonically increasing for o 0, 7 V0. 
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conditions satisfies the conditions of Case II. A necessary consequence of 


D) is se . 
saiesian ;(F"—@ 
| P21 G2 
P21 G” 
since J(A’) J(A”) 0, and the denominator does not vanish by 
hypothesis. Hence, 





ul I (F"/F’)_ ,  F" 

1—A’ F’24 G2 F’—G" — F?4G"? 

L+Ao (G"/@) 4 & 
F’24 GQ’ yg F’2+qg’2 


Finally, the relation J{(1+-A’)/(1—A’)! 0 leads to 


pr | ar j a” a” faa’ “hel \ 
Hi 7 Se ‘a 7 .) a : -)a( =. 3} J) 
f G I 21 (7 | A G wi 21 (}'2 


By a detailed discussion, to be published in full later in another con- 
nexion, it has been shown that the only solutions of (J) not already dis 
cussed which satisfy the conditions of the original problem are such that 
A’(u) constant.+ We must therefore have A(u) vu +, which leads to 
the yield function f(c) vo +-constant as in (I 1). However, the restricted 
choice of F and G@ in (11) led to only part of the complete family of 
trajectories for this family of yield functions. We see that the most 
general F and G satisfy the relation (a+1)F”+(a—1)G@” = 0. Hence, 
when a tI, 

F d* a, C1) G : w2+a,w—+c,. 
l+«a | x = di 
When a +1, we obtain exactly the family of nets for which the nets 
found by Carathéodory and Schmidt are the 45° trajectories. As we have 
already seen, the choice « = 0 falls under Case I. When c = 0, the nets 
are logarithmic spirals or their degenerate cases; otherwise, we have the 
nets given under (I 4). 

When «a takes on the special values | and — 1, the family of nets is seen 

to reduce to the degenerate subclass consisting of concentric circles and 


concurrent straight lines. 
6. Remarks on the linear yield condition 


We have found that out of the complete family of isometric nets 


belonging to the yield condition + vo +f, only the logarithmic spirals 


+ Further non-trivial solutions of (J) are given by the general solutions of the equations 
Fr" |F’ CF’+A and G’”/G@’ CG” — A; however, only the special cases which lead back 


to the above obtained solutions satisfy the conditions of the original problem. 


- 
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appear under (I1). In (12), on the other hand, we found the complete 
family excluding the logarithmic spirals, but only for the special yield 


equation 7 = constant. The complete family of nets for r = ao+8 corre- 
sponds to the only solution of (D) under Case IT. Whatever the value of « 
except for the special cases « +1, mentioned in the preceding section, 
7 = ao+f is consistent with the same family of isometric nets as the 
vield condition 7 = constant. 

The only linear yield equation not covered by the preceding is ¢ = con- 
stant. We may approach this case by introducing o = s(7) in section 3 
instead of r = f(c). In terms of this function, we obtain the simultaneous 
lene log V* log7+ F—G, (a) 

S(T) dr PF C (b) 


in place of equations (1) and (2). 

It is seen that the general mathematical structure of the problem in this 
version is the same as in the earlier one. However, in the present formula- 
tion, we can easily discover the highly exceptional nature of the ‘yield 
equation’ o = constant. In this case, and only in this one, both F and G 
ire found to be constant by virtue of equation (b), and hence (b) does not 
impose any restriction on the choice of the isometric trajectory pattern. 
From equation (a) follows 7 cV?, 

Thus, for this exceptional ‘yield condition’, we obtain a special case of 
iresult found by one of the authors many years ago in a parallel investiga- 
tion in the field of elasticity (3). It was shown there that for any isometric 
pattern a five-parameter family of stress systems can be found, consistent 
both with the equations of equilibrium and with Ao = 0, the equation of 
compatibility for plane elastic systems. We see that o = constant is a 
special solution of the compatibility equation; hence it is no yield condition 
in the proper sense of the word 

Finally, we can show that the linear yield condition is the most general 
me with which the Carathéodory—Schmidt nets are consistent. We see 
that A(F+G) must equal a function of ¢ alone plus a function of % alone 
n this case, and this is so only if A(w) vu+constant. The corresponding 


vield equation has been found to be 7 1o constant. 


7. Stress trajectory patterns consisting of congruent curves 
In this section we consider isometric nets such that each family of curves 
goes into itself either under a translation or under a rotation. It will be 


The cases «a l are indeed exceptions 


l, since the general system of equations (L,, Lg, 


nd the vield condition Ss paraboli ( i (a) a 
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seen that any one of these trajectory systems is consistent with the equa- 
tions of plastic equilibrium for a suitably chosen yield condition. 

In the case of translation let us assume, for simplicity, that its direction 
is parallel to the y-axis. Then obviously the direction angle # (measured 
from the positive x-direction) of either family is independent of y. Hence 
for such a family of curves, corresponding to the imaginary part of Z(z), 





, dZ 
o im | log | —cx+e, 
dz 
' dZ , ; 
and therefore, log — acz+-C. 
dz 


Here and in the following discussion capital letters denote complex con- 


stants, small letters denote real constants. If c 0, we obtain the trivial 
solution Z C,2+0,. (14) 
Ife +0, Z = C,e+,. (15) 


The corresponding isometric pattern is seen to consist of the curves 
cy = logecoscx+constant and cy = logsinca-+-constant. 

The case of rotation about the origin is reduced to the preceding case by 
the mapping z, = logz.t Rotation about the origin in the z-plane is seen 
to correspond to translation in the z,-plane parallel to the y,-axis ; hence 
the solutions in the second case are 


Z C, log z+C, (16) 
and Z = C, 2+. (17) 


The trivial solution (14) is, as we have seen under section 3, consistent 
with any yield condition. Family (16) is consistent with the general linear 
yield condition, as shown under section 5. As to (15), we find that it is 
consistent with the ‘hyperbolic’ yield condition mentioned under section 5. 
Finally, (17) is consistent with any yield equation obtained by setting 
a 1, B ci in equations (10) of section 5. 


8. Conclusion 

Apart from the exceptional case in which plasticity is replaced by 
elasticity, we find that a given yield condition is compatible only with 
certain definite families of isometric nets. In general, these nets are only 
the trivial Cartesian nets and nets consisting of concentric circles and 
concurrent straight lines; for three distinct families of yield functions f(0), 
however, depending on from two to four parameters, other farilies of nets 
are also found. These include family (10) (to which belongs the parabolic 
yield function of von Mises), a general linear function (thus answering our 


+ The authors are indebted to Mr. Daniel Shanks for suggesting this approach to case of 


rotation 
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oul question for each line segment of which the St.-Venant yield condition 
consists), and a general real power of a linear function. 
tion 
ired The last two yield functions can be obtained from (10) by limiting processes. 
nce Ak + KxT Dead 
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THE BENDING OF THIN ANISOTROPIC PLATES 
By EDITH ILLING (University College, Cardiff ) 
[Received 9 November 1950. Re-submitted 28 February 1951] 


SUMMARY 

A complex fourth-order differential equation is derived for the normal displace- 
ment of the middle surface of a thin orthotropic plate subjected to transverse 
pressure. The displacement equation is solved for the problem of a circular plate 
which is clamped around its edge and bent by a typical general force applied to one 
face. Particular cases of the general solution are used to investigate the bending 
of the circular plate under the action of (1) a uniform pressure and (2) a linearly 
varying pressure. The normal displacement of the middle surface is also determined 
when the applied force varies as the square of the distance from the centre of the 
plate. 


1. Introduction 

THE problem of the bending of a thin plate of isotropic material has been 
treated by several authors. Stevenson (1) gives detailed solutions for the 
circular, elliptic, and equilateral triangular plate under uniform pressure, 
and Seth (2) for the elliptic plate with a confocal hole. The general result 
for the deformation of the n-‘sided’ polygon is given by Stevenson (3), and 
solutions for plates bounded by cardioids, lemniscates, and other quartic 
curves have been derived by Sen (4). 


Few attempts have been made to solve similar problems for plates of 


anisotropic material. Morkovin (5) has given a complete solution for the 
bending of a finite elliptic plate of material having at each point a plane 
of elastic symmetry parallel to the neutral plane of the plate. Holgate (6) 
has investigated the problem of thin perforated aeolotropic plates bent by 
couples at infinity, and gives solutions for plates with (i) a circular hole 
and (ii) an elliptic hole, the material of the plate having two perpendicular 
directions of elastic symmetry parallel to the middle plane. 

Holgate’s general equation for the transverse displacement of an aeolo- 
tropic plate under normal pressure, which has also been derived by Huber 
(7), reduces to the equation given in this paper for plates of material with 
two directions of symmetry parallel to the middle plane. 

The normal displacement of the middle surface consists of a particular 
solution of the differential equation, with suitable complex potential 
functions added to ensure that the boundary conditions are satisfied. 
The conformal transformation used to solve the problem of the bent 
circular plate has two branch points inside the boundary of the plate, 


[Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 1 (1952)] 
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but this difficulty is overcome by constructing potential functions which 


remain uniform over the area of the plate. 


2. The displacement equation 

The plate is assumed to be of uniform thickness 2h, and rectangular 
axes of « and y are chosen in the middle plane of the plate, the axis of Z 
being at right angles to this plane. 

In all the problems considered in this paper it is assumed that the 
material of the plate has three mutually perpendicular directions of 


symmetry parallel to the axes of coordinates. For such material, 


= - if f, 
at $41 UL +819 YY + $43 ZZ 

Eo ey sz 
ae $41 TX A-Sy0 YY +853 ZZ l (2.1) 
Cc. y S66 Ly 


where the elastic constants s,. are the coefficients C../Il of Love's theory (8). 


With the approximate values of the strain components given by 


» 9 » 
Ow OW py OW 
_——, ( Z—., _- 2Z- , (2.2) 
Coz” coy~ CXLCY 
inserted in (2.1), integrating, we derive: 
or h 
2h3 o*w ' ae ae 
811, 4 +832G,+-83 | ZZZ dZ, (2.3) 
5 Ca” . 
: h 
h 
2h ow Y y 7 l, wT, ’, 
= $10 G1 +80 G.+8., | ZZZ dZ, (2.4) 
» oY 7 Se oiler 
Y h 
1),3 O-u 
Seq fi, (2.5) 


» . 
» CLO" 


where the stress couples G,, G,, H, are defined in the usual way. 


h =s 
Solving equations (2.3)-(2.5) for G,, G,. H,, neglecting | ZZZdZ in 
I, 
comparison with G,, G,, H,, we find 
; 2h?(  d*w cw) 7 
Gy Ste oe Te or (=.6) 
3 cx? cy” | 
2h? ( d*w c*uw| ™ 
Gy a ite 5 1 22 o(? (2.7) 
3 I Cx? cy") 
+h* o*u 
i, H, ee (2.8) 
2 a 66 a 
o CXCY 


W here 
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From equations (2.6)—(2.8) expressed in terms of the complex variable 
z= x+y, and its conjugate 7 = x—iy, we derive the complex stress 
functions 


; 2h? ( 02w C2w 
G,— G,— 21H, aaa (or + C99 — 2C19—4C gg) —j + 2(€43—C 2) —=+ 
3 | 02" 0202 


O2w 
+ (€yy +-€g9— 2€y9+ 408) <a} (2.10) 
O52 


(2.11) 


The essential equations of equilibrium for the plate bent by transverse 
pressure F' per unit area are 


oH, eG, 


mone: AD 8 0, (2.12) 
Ox cy 2 
CG. of ; 
lone a on @, (2.13) 
Cx oy 
C N, 4 ON F 0. (2.14) 
Ox oy 


From equations (2.12) and (2.13) written in complex notation we obtain 


N,+iN, = —(G,—G@,—2iH,)+—< 


Cx Cc 


(G,+G,), (2.15) 
and, in complex notation, equation (2.14) becomes 

(N,+72N,) + —(N,—1N,)+F = 0. (2.16) 
Substituting in equation (2.16) for N,+7N, and its conjugate N,—iN, in 
terms of the stress functions G,—G,—2iH,, G,—G,+2iH, and G,+G, as 


given by equations (2.10) and (2.11), we derive the differential equation 


for the normal displacement w of the middle plane in the form 





‘Aw = Aw Cfw Clw : Cw 3F ee 
X4 {| —__ 1 Any. ae t+ | + 203-5 -- = = 9, (Z.8¢) 
oz* — az! “\oz80z  020z 0270z2 =2h3 
where a om 
XY C11 T&22— “C42 CE6> 
Xo C11 — a2 
> Qn Ip > 
X3 001177 3c 227 <"y277 AC 6g. 


This equation can be obtained from Holgate’s general displacement 


equation by taking both c,, and c,, to be zero. 
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3. Solution of the displacement equation 

The solution for w consists of two parts. One is a particular integral of 
the differential equation with F appropriately interpreted as a function 
of z and 2: this we denote by Q,. The complementary function is of the 


characteristic form 


1 

u > 2.4.) (3.1) 
m1 

where z,, = 2+-A,,Z, and A,, As, Ag, Ay are the four roots of the reciprocal 
acai (A441) + 4ev(A3-+-A)-+ 2avg AZ = 0. (3.2) 
Since two of the roots are reciprocals of the other two, we may write 
h, = L/A,, A, | /A,, where A,, A, are real or are complex conjugates, and 
ire chosen to be the roots of modulus less than unity, i.e. |A,| < 1,|A,| < 1 


The subsequent analysis is simplified by taking A,, A, to be real. The 
normal displacement of the middle surface then takes the form 


0(2,2)+ ¥ {Qn(%m)+Qin+2(Fn/Am)}- (3.3) 


) = 0,(2,/A,) = 2,(2,). (3.4) 
und, similarly (2 ,(Z5/As) (2.2). (3.5) 


Thus the expression for the normal displacement reduces to 


w = 0,(2,Z)+ > {0,,(%n)+Qn(Zmn)}- (3.6) 
m= 

The assumption that A, and A, are real does not hold for many aeolo- 
tropic materials, but the method used for real roots can easily be extended 
to the case where A 


Ne gf 


;, A, are complex conjugates chosen so that |A,! < 1, 


1. The stress functions in terms of the displacement 
The complex stress functions can be expressed in terms of the displace- 
ment as follows, dashes denoting differentiation with respect to the 


irgument of the function 
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G,+4, {(1 +-A?)(QY +O) — (1 +.A2)(Q3 +03)! 
6A, A. 
L6h3c " : i. 
of fA (Q)4+-Q0)) A,(Q3+Q3)} 
Zhe, (eFQ eQ) e7Q) 
2 ol ‘’ =) (%3— Oy C66 4 > (4.2) 
a 62? ez? 20z | 
N,-+-1N, 


h3x.,(Ay—A,)(1—Ay Av) pc aaiee ia pS tas ei 
Tals AEA 9) 4 92)(,, 2" —O") — (1-8) (Ay Q"” —O 


6A, A 2))— 
= 
D3 3 3 a3 3 
2h f OQ ee: Q2, O°Q6 CQ | 4: 
ee Sy ae Stig a pe ge) «= C2) 
3 |" é 02702 0202" Oz | 


The normal shear, reckoned across the whole thickness of the plate, for 


an are element ds of any curve drawn on the middle plane of the plate is 


_ dy ax 
lids “ds (4 4) 
which, in complex form, is 
ee Te <4 serge 
ge te a 
The couple acting on the element ds is 
, ., dx : ., dy 
G ; G : 
(G,i H,j) — (H,i Ie 
ij G.—iH dz | = iH. 4i@ e =) 
ai, day NG, del 
i | ; , : dz ; ; dz) : 
(G G.— Riz G,+G, aa (4.6) 
gf Mag, = 


where i, j are unit vectors in the directions of the (2, y)-axes. The tangential 
component of the couple is 


: dz\* ; dz\? dz dz 
LG, —G,+2iH) } WG —G,—2iH ( va uggs &. am 
. ; as WG. . ss ds 2 . FP ds = 
and the normal component is 
i aa : a 8... ; . (dz\ ’ 
riba Gy 21H} d rise (ry 2iH}( ; (4.5) 


8 ds 


5. The characteristic stress potentials 


The normal displacement w of the middle plane has been found to be 
of the form 


9 


uw 42,(2,z > 492, (Z,,) +2, (Zn) ts (5.1) 
m=1 
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where Q, is a particular solution of the quartic (2.17) in which F is ex- 
pressed as a function of z and Z. The problem is to derive suitable functions 
Q. (2) to ensure that the boundary conditions are satisfied. The method 
of solution used is essentially that described by Livens and Morris (9). 
The displacement equation is solved for a circular plate which is clamped 
>) round its edge and bent by various transverse forces. The conditions 


which must be satisfied at the clamped edge are 


ow Cw 
we (0. 0 ( ‘ 0) (5.2) 
C2 Cz 


[he boundary conditions have been stated in this form by Stevenson (3). 


3) We take as the boundary of the plate the curve 7 0 in the net defined 
™ (le C é in. (5.3) 
Is Thus. on the boundary 
t) r,,z = a(t+A,,/t), (5.4) 
wtnien 2 exp! 
We wish to determine f,, as a function of z,, which reduces to 
5) } 


expl 1€) 


the boundary. If Q,, is chosen as a function of ¢,,, both Q, and Q,, can 

ve expressed in powers of ¢on the boundary, and comparison of like powers 

in the boundary condition equations will determine the unknown 
efficients in the functions Q2) (7 


Equation (5.4) has two roots 


fe 4 i= r,, , (5.5) 
6) 2a A 4a 


} one of which reduces to exp(—7é) on the boundary. The roots become 

- jual at the points z,,/2a vA,, Which, if A,, 1, are both inside the 

te. Thus there are two branch points of the transformation in the area 

=) nder consideration, and the roots ¢,,,, t,,,. given by (5.5), are non-uniform. 

. One branch maps the interior of the plate in the z-plane on to the ring 

| ice | . \A,, in the ¢-plane, and the other branch maps the 
terior of the plate into vA ae 

.5) \lthough the roots ¢,,, ¢,,. are unsuitable in themselves as potential 

inctions, t” ” is uniform over the area of the plate since, from (5.5), it 

msists only of a finite power series in z,,. Since ¢,,,t,,, = A, the uniform 


be potential funetion can be written in the form ¢",-+-A%,/t,, where t,, denotes 
he root of (5.4) which reduces to exp(—i&) on the boundary. 
. Thus suitable functions Q, (f A" t) have to be added to Q, to 


) 


( 
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determine the complete displacement of the plate with the prescribed 
boundary conditions. 


Note that ét,, l ét,, G. (5.6) 
éz A, ez  a(f,—A,)’ 7 
l ¢ (. c a i 5.7 
A, © cz a(t? —X,,) _ 


and, in general, 


(0 | “ l ¢ (0 | in) MY in 1_| x tn 3 ee Nn : | rn 


m jn a\ m m*m yn 3 yn 1 


m m m 


Cn 


(5.9) 
6. A general solution for the clamped circular plate 


Consider the circular plate defined as in (5.3) bent by the transverse 


force >] 3 ' 
2h°a,(n—1)! ,,, > ; 
h oad — {K, 2"-4*+-K,, 2". (6.1) 
3a"(n—A4)! 
With this value of F inserted in the displacement equation of the plate 
Cw cw Cw cw > eMw 3F rv 
| 4 =4 1 tx, 305! aah) ““$7 sense" h3 V, (6.2 
| OZ Cz O2°Cz C202" 02°02" 2h 


we obtain the particular integral 
{,, (K,, 2"+K, 2")/na". (6.3) 
The boundary conditions for the clamped plate are 


. 7 cw 
(i) ow 0: (11) 0. (6.4) 


Cc 


Assume the complete normal displacement to be 


w (K 2®+kK =") /nan \ % m,n—2) (" 2r_) “m 
" Dap 


y A ‘ m op 
Lag, Valent |, 29 — Qs o> 
1 


r=0 
- n—2r 
Amn 2r pr 2r “\m | (6 5) 
Dar m 7; Dp 3 ” 
n— 2? lias 


the upper limit for r being (n—1)/2 or (n—2)/2 according as n is odd or 


m 


even. 
Then, 
9 
. -_ 1 cn if “ f n— 2r—1) 
ou K eg l \ IX a jn—2r—-1 \ fn—2r-3 “m ly 
~' / ) mn—2r m “Sm “m T eee] n > 1 
Cz an a Lew | Ly otis 
m=1 ‘r=0 . : 
9 
|< “ ~ a : \n— 2r—1) nl 
S | he Ie- emt) Pa-e-8s. 4 \_ (6.6) 
Jf m myn a m mom n dr l | 
a peal \ — wk a 
m r=() 
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On the boundary, ¢t, = t, = ty? i, 1 = t = exp(—7€), and to satisfy the 
boundary condition w = 0 we must have 


yp TAR Gyn +O n+ABG,,+K, = 9, (6.7) 
Bs n—er TAT 7G, nor +e n—op tAZ Gene = 90 (9 1, 2,...). (6.8) 
Since ew/éz 0 on the boundary 
A, »+At a, Ay ,+Az a, ,+K,, 0, (6.7) 
Y Mu Ginn tM Ginn tA Cinn—2+ Am” inno 
d 1 
ee ee 0 (r Pe ene (6.9) 
With r 1, equations (6.8) and (6.9) become: 
Oy ng AT Ay ng +o, n-2 AQ "Ge n-2 = 9; (6.10) 
Ay 44, +At ay Ag by +A gn +41, n-24 
NG, n—-2 +e n-2tAZ Ae n-2 0. (6.11) 
Hence A, 44 ntl 4, Ay As » +-Ap ld... 0. (6.12) 


Solving (6.7) and (6.12) and their complex conjugates, 
a, a, As, ay (K,, T K,,) . 
A,(1-++AR-?) Ay(1-+-AR-2) (1 + AR-2)(1+-A8)—A, (1-+AR-2)(1+.08)’ 
(6.13) 


and 
a, a, ay 7 (K,, —K,,) 
\,(1—Az-2) A(1—AR-2) —-Ag( 1 —AR-2)(1 — AP) —A, (1 —AR-2)( 1 — AB) 
(6.14) 
With 7 2. (6.8) becomes 
a P rt Ma, st+-Gen-~ Ay “... PF 0. (6.15) 
With r 2 in (6.9), and (6.15) taken into account, we obtain 


2 n—27 2 n—27 n—-3q 
Apa, Atay NB 9,» +AQ "Ge n+ Ay Gy ng FAT Gy» -24 


Ag@. ,-0+Ag-G,,,-2 = 0. (6.16) 


Solving equations (6.10) and (6.16), using the values of the coefficients 


a,..a,, given by (6.13) and (6.14) 


(il, 9 a, 2 as . a . A, AK, 


K,,) 


’ 


1+-\n-2 12-2 Ao( 1 +-AZ-*)(1 +-AP) —A, (1+-AP)(14 Ny)’ 
(6.17) 
ind 
a a A, A(K,,—K,,) 
1—)z-2 ie Ag( 1 —AR-2)(1 —A”)—A,(1—AR-?)(1—AB) 
(6.18) 
Thus a \ 


a2 A, 413 fe n-2 Ay Ms ,; (6.19) 
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and it can readily be seen from equations (6.8) and (6.9) that the remaining 
coefficients @, ,,-2-, 42. n—2, (r > 1) are all zero. 


Thus the complete normal displacement of the middle surface is 


P. coe -—. { n 7 f n 
* K,,2°+K,,2" . ay mn P\ . ai. zn rN 
+ 1 2 iT Sr 

na” n tt n tt 

n 7 n n—2 

2.n yn AS As» jn 3 Ay 41.» jr-2 Ay 
eae 1 ‘ty ie si ee 1 —— 
n im) "nm iz} n—2 in? 


4 n—2 -_. 4 n—2 
A, ual te 2 ri “a Sao (eg 2 AS As Gs,» jr 2 ro 
ap thas a 9 \"2 (ia 312 | Sas)? 
n—2 5" n—2\ ; n—2 ts~* 


where the coefficients a, ,,, @,, are given by (6.13) and (6.14). 


It can be seen from equation (6.14) that the coefficients a, ,,. a,,, are 
real when K,, is real. The complete displacement is then 


r (on _| =n , nn = An 
a. K,,(z ) \ a - ti - 4 


7 ” ' n 
na n t ty 
n n f 9 n—Z 
Aonfin AB, in, %B An [yn 2 ri wee, Ms 
ty Ton ty inf ) fy T on—2 | fi T “in—2 
n ty ty n—2 i“ tithe 
Ao As «ae .. ae 
sel = ay = (6.21) 
n ») - {” 2 - fr 2 
-_ » » 


where 
A.(1-+Az-2)K, 
in No( 1-EAR-2)(1-+AP)—A, (1+ AR -2 
a A, 14 vA; 2\K, 
om Xo(1-+EAR-2)(1-LAP)—A, (1 +-AR-2)(1-+A2) 


rN) 

(6.23) 
Particular cases of this general solution given by equation (6.21) are now 
considered in detail. 


7. Plate bent by a uniform pressure 
With » 


uniform pressure 


4 and K,, real, the transverse force defined by (6.1) is the 


Sh3.x 
1 K - 

\ ). (7.1) 

ai 1 / 

The normal displacement of the clamped circular plate corresponding to 


this applied pressure, found by putting » 4 in equation (6.21), is 


*) (a r; L a+F) 
= 1 = 
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y where , 9 ms 
? Bs K,(1+-A)A,/Q, (7.3) 
dy, = K,(1+22),/Q, (7.4) 
with Q = X,(1+A3)(1+-At) —A,(1+-A2)(1 +09). (7.5) 
Since t Am/t, (2+-A,, 2)/a b+ Are/ Em (z+-A,,,z)/a, (7.6) 
me (1+-A?, )(z?-+-27)+-4A,,, 22 — 
we have i + t° = Am) = yr *\n 4X,,,; (7.3) 
E i a* 
and 
1 1 
/4 Xm 74 \ 
1 74 
)) 
4 ae 2 Viste 3 2 ~252) 
(1+-Aj,)(24+-2 4A, (1 +AZ, )(23Z-22' 12X°, 22z?! — 
€ at 
7 i 9 9 9 =) 9 ~ Oo 
((1+-A5,,)(27+-27)+-4A,, 27'+4A2,, (7.8) 
a 
Substituting these values in (7.2), with the expressions for a, 4, a, 4 given 
by (7.3) and (7.4), we derive the displacement in the form 
. 3A, A. K, 222z2 22% 
L+- 2A, A,+-AZ AZ +-AZ+-A3\ at a? 
l) : 
3K, Xy a2 y2)2 
2a4a. 
) - 
op at 
(q- a (4 g 
L6hF x. 
3) Shi, | 
where ‘44 and 0 r—ae<a. 
as 
W 
For isotropic material, 
, 2u(A+p)/(AS 2), (7.10) 
and (7.9) reduces to 
he , \ 
3 1, Zu > 9\9 ) > \9 - 
u : “i (a*—r*) f (a*—r?)?, (7.11) 
512uh3(A+-p) 64D 
where D 3uh3(A+p)/(A+ 2p). 
to 
This is Love's result for the problem for the isotropic plate (10). 

It is interesting to notice that the anisotropy of the material has no 
| effect on the general shape of the plate, the law of variation of w being 
| of exactly the same form as for the isotropic plate. 

The deflexion w takes the maximum value 3pa‘/16h3x, at the centre of 


t] 





ie plate and decreases to zero on the boundary. 
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The tangential bending moment couple on the element dé of the circle 
‘- 
Z re~*s 18 


: a de\2 — di\2 , , dz dz 
1(G, G, + zim) (Te) + L(G; G, zit) (7; 3(G, + Oe) dé’ 


Using the expressions for the complex stress functions (G,—G,+ 27H,), 
(G,—G,—21H,), and (G,+G,) in terms of the displacement given in 
equations (2.10) and (2.11), with w given by (7.9), the tangential couple 
on dé is found to be 


P fa,(z4-+24) + 20(322—a2)(22-+ 22) 


16g. 
€ ae ae Wos . ne oe 
2 (a + 8C¢4)(22 —@*)2z + 2a9( 222 —a*)z2} 
pr 
{x, rcos 4€ +- 2a,(37?—a)cos 2E + (a+ 8¢44)(a2—1?) + x4(2r?—a?)}. 
SX 
(7.13) 
On the edge of the plate, where 7 = a, the couple has the value 
par et a} - 
— ja, cos 4€-+- 4a, cos 2€+- ag}. (7.14) 
Sag . 
The stationary points occur where sin 2é = 0 or cos 2& y/o. At 


points on the boundary where sin 2 = 0, the corresponding values of the 
tangential couple are 


., par ae ee 
(1) —}ay+4ao+ ag; (E 0, 7), (7.15) 
Sag 
- a4 ‘ wv 3 rf 
(11) l 1X, —4a-+ a5} (¢ ; } (7.16) 
Sag ‘ : 2 2 
At the point on the edge of the plate where cos 2é v,/a,, the value of 
I Lu I 2/%y 


the tangential couple on the element dé is 


14 202 ee 
- (» x4 — 3), (vie 
Og Xy 


The normal torsional couple on the element dé is 


a ’ ’ aa dz\? d , Y 2; dz\* 
ih Gs | 2 (5 - rides Gs, >it) (54) 


Pr sin 2¢{a, r2cos 2€+a(2r2—a?)}, (7.18) 
X3 
which vanishes at the stationary points of the tangential bending couple. 
On the boundary, the normal couple takes the value 
pat , 
= ( 
405 


Yo +a, cos 2&}sin 2€, (7.19) 
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THE BENDING OF THIN ANISOTROPIC PLATES 
le The normal shear on the element dé of the circle z re-€ is 
Pr Sree pr? 
(N,-+-+N,) — (N, —2N,) = {3a, cos 2€+a,'. (7.20) 
* ~ @e 2 “dé as ' " , - 
The normal shear is zero at the centre and increases to 
9 
“) pa- ’ » _ 
\ {3a, Cos 2€+- ag}, (7.21) 
1) 23 “ 
i ; ; , 
; at the edge of the plate. It is a maximum at the edge of the plate for 
le Ay! a, F 
& = 0, 7, and a minimum for & 7/2, 32/2. 
8. Plate bent by a uniformly varying pressure 
When n = 5, and K, is real, the transverse applied force defined in 
(6.1) is a uniformly varying pressure: 
F 16h3x, K;(z+2Z)/a> = pox where py = 32h*a, K,/a°. (8.1) 
)}. i ' ; , 
rhe corresponding displacement is 
3 ? 5) 3 3 
re K; (25-4. 55) 4 “(4 A\ 75 Ai AAs 34.01 4 F384 Ay 4 
5a ‘i. 2° = Boe bia 
5 5) 3 3 
1) Qe5[5_) 42 75 Ay Ay My5 pp A3 73 2 (8.2 
a "eS Bit tee 
At where ins K,(1-+23)A,/R, (8.3) 
he ay, — K,(1+23)A,/R, (8.4) 
R = X,(1-+A3)(1+-A2)—A, (1 +A9)(1 +23). (8.5) 
5) Also 
a, 
= + bin + 5 
6) tm f; 
l S - ~ ¢ ‘ -* - 9 - 9-9 
ad atl AS, )(z> +25) + 5A, (1 +5, )(z3+-23)zZ + 1OA2,(1-++A,,,)(2-+2)2727} 
5A , ‘ = : aa 5d? . 
(1) AB )(z3+-23) + 3A,,(1+A,,)(2+2)2z} + —™ (1+-A,,)(2+2) 
7) a a 
und (8.6) 
\3 x \3 
{3 | m | 3 mm 
n 73 ” #3 
l . ire —-_ ™ ana - 
311 A®_)(z2+-23) + 3A,,(1+A,, )(2+2)2z} - (1+-A,,)(2+2). (8.7) 
a a 
18) eee ; 
Substituting these values in (8.2) 
le 2K, A, A(z +2)(27Z2— 2a*2zZ+-a4) 
uw 
| a®{2Q, Ag+ 1—Ay+-AZ—Ag—A¥ Ag+AZ—A, AZ+-AF AG} 
19) | re ie (5.8) 


16h3( X3 
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For isotropic material, 
).) 
on u(A rp 
9 were = ®, (8.9) 
A+-2p “ 


and (8.8) reduces to 


(A T 2) po X 9 9\9 Pot 9 9\9 
v | (a2—r2)2 (a2—r2)2, (8.10 
' 512wh3(A 1) 192D ) 
where D $uh3(A + 2) (A 2). 


This is the result given by Love for the isotropic plate problem. 

Again, the law of variation of w is of the same form as for the isotropic 
plate. The deflexion is zero at the centre and at all points on the y-axis, 
and is a maximum at the points z La/wv5. 

The following results are determined as in the case of the plate subjected 
to a uniform pressure: 


(i) The tangential bending moment couple on element dé of the curve 


" re—€ is 
if ‘ € -\9 ) -9 > =-9 os 
ro (a4| (322 OO" )\S—Z)° 4-2" -- 2° 2-4-2" 2z)|-4 
£8(a, 209) 
2u,(322 —2a?)z7+ 2a, (322 —2a?)(z?+- 22+ 22) + 22(22+4 22—22)] 
1 6C¢6(2 a*)zz' 
yA ) 9 9 " 
— {a4[2r? cos?2é-+ 2(r2—a?)cos 2€— 37? + 2a2] 
» 
24 (a, 2x5) 


¥3(372— 2a?) Zao) 2(2r? a?)cos 2€-+-r2 a*| 8Cg,(77—a?)}. (8.11) 

The couple vanishes at all points along the y-axis; on the edge of the 
plate, it takes the form 

Pp @°COs & 


24(a3-+ 2a) 


‘ 
( 


y, cos 4€+- 4a, cos 2€+ az}. (8.12 


(ii) The normal couple on dé is 


p,(2—Z) _ . 9 e) a6 ~ ; > - 9\.5 
ro {(2+2)?| a,(2227—2a?+- 224-2?) +4 2a,(32Z7— 2a?) |+ 16c,,(22—a?)z7} 
96(a, 2a) = 
9 
gy hY o¢ ) o¢ 9 oo 9 9 , ‘ 
- Fo" {cos*€| «, (27? cos?é —a?) + «,(3r?— 2a?) |+ 2cg4(r?—a?)}. (8.13) 
6(a3-+ 2a) : 
It vanishes at z Lia and at all points along the x-axis. The normal 


couple on an are element of the boundary curve of the plate has the value 


», asin E cos7é . 
Po - = 1a, COS 2ZE+- ao}. (8.14) 
2M) “ 


(iii) The normal shear on the element dé is 
Po(2+2) ‘ 9 9 


1 30,(22—2" 
48(a,-+ 2a) ' ” 


)— 3a.(422—2a?+ 322+ 32?) —a,(927—2a*)}. 


(8.15) 


A | EE LE OO 








valu 


9. ] 


The 


equ 


wh 
wit 


ant 


Su 


in 


di: 
1( 


m 











9) 


the 


12) 


nal 


ue 





Nee ores - 


~~ 


ro - 





| 
| 
| 





THE BENDING OF THIN ANISOTROPIC PLATES 


25 


[t vanishes at all points along the y-axis, and on the boundary has the 


value 

1, A°COS & ' “ i , ™ = . . 
i : : 1 3X, | | 2 cos 2€) 6a,(1+-3 cos 2€) i Net. (8.16) 
24(ae+ 2Za,) . 


9. Plate bent by a quadratic pressure distribution 


When n 6, the applied pressure defined in (6.1) takes the form: 
; 10K. h3a, SOK, h? , 
I Se 7% (224 52 to 0 M1 (72 2). (9.1) 
a® a® : 
The corresponding normal displacement of the plate, derived from 


equation (6.21), is 


r (B =f { 6 \ 4 4 
ae Ko(2°+-2°) “slit a, e, | “Sule | A; 1 f44 Aj 1 
ba' os = oe .. . > ae 
G96] 46 A3 76 | AY Az Mo 6 4 AS #4 A} 9° 
ls : te =< fn ts =|, (9.5) 
Rieke Ses : i" €°”’a 
where a A,(1+-AZ) K,/S: Ay = A,(1+A4)K,/S, (9.3) 
with S = A,(1-+-A$)(1+-A8)—A,(1+Af)(1 +8), (9.4) 
and 
X' ”? : = ee 
“pe : c —- f(1+A® )(z8+-28) + 6A,,(1-+A4,)(24+-24)22 
f' e a° 
1 5A2, (1 +-A2,)(z2-+- 22)z2z2+4 4008, 23Z3) 
6A - 9 9 -9 > 9-9 
m {(1-+-A5,)(z4+-24) + 4A,,(1 +-Aj,)(2?+-2?)22 + 12A,, 2727} 4 
a 
9X2 woe . 3 - 
mf(1+-A2,)(z?+-27)+ 4A,, 2Z!—4A5,. (9.5) 
rt Aa 


Substituting from (7.8) and (9.3)-(9.5) in (9.2), the displacement is given 
in the form 
LOAZ AS K, 
30%} 2A5 AS gt (1 AF) ( | dS) X45 


{3770 COS 2é 2(47r° A*)x9}(r* a*), 


(9.6) 


It should be noted that the anisotropy does not affect the form of 


displacement of the plate until the applied pressures become quadratic.t 


10. Pressure proportional to the square of the distance from the 
centre of the plate 
The deformation of the clamped circular plate is now determined for 


the applied pressure FR pee. (10.1) 


Solutions can be obtained in this way for higher powers of n, but the analysis becomes 


more complicated as n increases. 
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The displacement equation takes the form 
py ok Lp fee Hen Bus 
Ctw = fu ofu O*w \ . ofu 3p2z 
a actos BENS Aaa aciel F a oo 
oz* oz4 “\ 02302 e a 2 


of which a particular solution is 


Pp Pore si ee F 

Q 3(2°z-+-22z° “(28+ 28)}, 10:3 
. F80/%, | ( ) : ( } ( 3) 

(i) Consider first Q, PM2 (61 36), (10.4) 


240h3 x2 
Quoting (6.21) with K, pa®a,/40h3a5, the displacement corresponding 
to this part of the particular solution is found to be 

Py ri AZ 


w 39 EY {3r7a, cos 26 — 2(4r?—a?)a,'(a?—r?)?2, 
L2h3a2 {207 AR agt(1+AP(1+A$)ay}o * , ” 
_ : (10.5) 
(ii) The remainder of the particular integral, 
y oo _ 
Q, f (2°z-+-22°), (10.6) 


160A3 x, 
is not one of the standard functions, and the complete displacement 
corresponding to it must be determined from first principles. 
The complete displacement must satisfy the boundary conditions 


Cw 


(i) w= 0, (1i) 0, (10.7) 
C2 
at the clamped edge. 
Write f 
aa 6 2 4 ’ 2 
ue K ( gs 22°) e 7 | A m (« a Xin | B,, (0 * Xn rl ( m (x a A | 
a, 6 . On 4 a. e. 2 = G. 
2. A m7 Ni B Z Xt. CA 7 te 
( = (7, | mi | rm <i | x I a e ! — : ( | 0.8) 
1 | 6 ip 4 vin - bin 
) 
where a —. 
160/3x, 
Then 
: 2. 3 5 
Be — KB + 2) + > (tt Att Ml EL 
- m=? . bm tn m 
2 2 3 ’ 
~ { Bi, f | A t eo Xn | Xn | ( mit i An | - 
a | a m m*m ‘ 13 a m t,, | 
~ {Am A m 75 r 73 2 i - Nn r>,, | 
a \ a m mom mm ‘. nd ad | 
2. Df. 2 3 ‘ta 
S (An B,, nd N Asn a a Xn ! Xn | Xn m ‘. aks Aw | ; ( 10.9) 
= \ a bm bmn a . ty, 
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Since A is real, the coefficients A,, A,, B,, B,, C,, C, are also real. Thus, 


the boundary conditions require that 


(1+-A})A,+(1+A8)A, = 0, (10.10) 
(1+-Aj)B,+(1+A3)B,+4Ka* — 0, (10.11) 
(1+-AV)C,+(1+A3)C, = 0, (10.12) 


A,(1+-Aj)A,+A,(1+A})A,+ (1 AN) B, +-(1- AS) B. -5kas 0, (10.13) 


(1+-AZ) B,+A,(1+-A3) B44 
+-(1+-A¥)C, + (1+AZ)C, = 0. (10.14) 


Equations (10.10)—(10.14) determine the coefficients A,, A,, B,, B,, C,, Cs 


uniquely as 


1, (1+A§)Ka®/S; A, (1+-A$) Ka®/S, (10.15) 
B f P(1+-A$)— 4A,(1-+.A2)| Ka®/¢ 
2 o( 2) "i é \ (10.16) 
B, {P(1-+-At)—4A,(1+.22)} Ka8/Q J 
( 5A, A,(1-+-A?) Ka®/¢ 
tig Pe Y \ (10.17) 
c. 5A, A.(1-+-A?)Ka8/Q J 


where 


S = d,(1+-A$)(1+A8)—A,(1+A9)(1 4-28), 
P (A? — A3)( 1 —A? AZ)(1+-AZ)(1+-A8)/S, 
Q = A,(1-+-A3)(1+A#)—A,(1+-A2)(1 +29). 


Inserting these values for the coefficients in (10.8), and expressing 


‘s A > 
f? my fe . (m 1.2 a 2, 4, 6) 
e = 


in terms of z and 2 as in previous examples, we obtain 
K «,(r?—a?) ; 
2A? AZ xg t-(1-+-AP)(1+A8)a, 


(8A; 


AS (ag— a, ) Xs Xo ° 9 Df of ‘ 
pe SN * | 2a? 2 (3r2—a?) + 5r?(r2—a?)cos 2€] — 


2 (AZ +A, Ay +-AZ)(1 +A, Ap +AZAZ)(10r4 — 8a?r? 4 a‘)! 


Ka 9/..9 9 9 -_.9 
A ! q?(r2—a?)(6r2—7a?). (10.18) 
X3 
The complete normal displacement of the plate produced by the force pzz 
is then obtained by adding (10.5) and (10.18). 
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Thus, we derive the final solution 


d2 A2(r2—a2) p 


w 
240h3.x,{2A2 AZ ag+ (1-+A4)(1-+A$) ay} 
j 9.9 4 ) ~ 9 9 9 1 2a 9/9,.9 9 | 
x 1(109r S8a*r-+- a*)( Zag x1) 12a,| 5r*(r-—a“)cos ze -“a*(3r-—a*) | 
; ¥3 
"a a 
. (6r-— 7a*)(r-—a“). (10.19) 
160h3 x, 
For isotropic material, the displacement reduces to 
pa (6r2— 7a2)(r2 2) 10.2 
w yr-— (a-)(r-—a*), (10.20) 
1920D 
where D Suh3(A+-p)/(A+ 2p). 
Conclusion 


Solutions for some of the bending problems of the clamped circular plate 
have been derived by the method of complex potentials, despite the 
existence of branch points in the area under consideration. The results 
indicate that the anisotropy of the material does not become evident 
unless the distribution of applied pressures is of quadratic or higher order. 
In the latter cases the analysis is very involved, but unique solutions for 
the deformation of the plate can be obtained. 

Finally I should like to acknowledge with gratitude my indebtedness to 
the late Professor G. H. Livens, who suggested the problem and subse- 


quently gave me much helpful advice and encouragement. 
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ON THE COMBINED TORSION AND TENSION OF A 


PARTLY PLASTIC CIRCULAR CYLINDER 
By F. A. GAYDON 

(Dept. of Theoretical Mechanics, The University, Bristol) 
Received 17 July 1951] 


SUMMARY 


Particular combinations of twist and extension of a solid circular cylinder are 


nsidered. Tl \euss equations are used throughout and these are integrated, for 
erent cases, to give the shear ss and tension in the plastic material. It is shown 

it the stresses rapidly approach their asymptotic values and are within | or 2 per 

, f these values when the plastic strains are still of the same order as the elastic 
trains. A more general case, in which the twist and extension are such as to make 
ratio load to torque constant, is solved numerically. Finally the residual stresses 
valuated, after partial unloading, for a bar which has been twisted and extended 


mstant rat 


1. Introduction 


[HROUGHOUT this paper the Lode parameters p, v are assumed equal. In 


this case it has been remarked by Hill (1, p. 75) that if, in the combined 


torsion and tension of a fully plastic cylindrical bar, the elastic component 


f the strain is assumed negligible, then the only non-zero stresses are o, = o 
md + 7. and there is ho \ ping of the cTOss section. The velocity 
esolutes u,v. win cylindrical coordinates are 
ly Orz IE 
9 ie i 
here the origin of coordinates is at the centre of one end. which is fixed. 
d the z-axis is parallel to the generators of the cylinder, / is the length 
the bar, 6 the total twist. and a dot denotes the time derivative. 
We now show that the above velocity resolutes may be used when the 
is only partially plastic provided Poisson’s ratio is taken to be }. 
\ssuming that there is no work-hardening, the Reuss equations (1, p. 39) 
r the plastic region are 


I/ a dx l 


/ l do. (1) 
»/ ) 6G 
’c AX l 
le do. ( >) 
2 3G 
a | 

dy, r dx : do, (3) 

d> ly 0. 


[Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 1 (1952)] 
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where dA is a positive quantity varying in space and time. The von Mises 
criterion of yielding is 


o? +37? Y2, (4) 


where Y is the yield stress in tension. 


2. Combined torsion and tension with constant relative extension 
and twisty 
Let A denote the initial length of the rod and a the external radius. As 
the load and torque are increased from zero the stresses o and 7 increase 
until at some moment of time equation (4) is satisfied locally. While the bar 
is elastic, o is independent of r and 7 is directly proportional to r. It there- 
fore follows that the material at radius a becomes plastic first. Subse- 


quently, when the inner radius of the plastic region is c, the stresses in the 


elastic zone are ss 3G Ini(I h)) . | : 
r Ci. b 
T Gré/l | (9) 
Co 31 
and we note that pint! h). (6) 
cz r 


Eliminating dA from equations (2) and (3), and using (4), we find 


3 dl o YY?! 1dr 


yr dé T G ror dé 4) 
for ¢ ? a. 
If the relative rate of extension to rate of twist is constant. then 
3dl/dée an. 
say, is a constant and equation (7) becomes 
an o Y*l 1dr 
(5) 


r T G ror dé 
If the value of 7 in an element is constant, o is constant from equation 
(4), and a solution of equation (8) is 


(9) 


G ——- T ae (Cr <a). (10) 
(a°+ or a*) (a* 1. 476 a-) 


Equations (9) and (5) are compatible at the plastic-elastic boundary since 
it follows from equation (9) and the definition of « that 


C ax 3/1 3 
we ( “) ~ ( jin h) 
T r r Q rd 


to the usual approximation in the linear theory of elasticity; o/z is therefore 


+ This solution was communicated to the author by R. Hill. 
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continuous at the plastic-elastic boundary. This, together with the fact 
that equation (4) is satisfied at both sides of the boundary, proves that o 
and 7 are separately continuous at the boundary. 


The load L and torque 7’ needed to maintain the constant relative rate 


3c°\1 
(a?-+- 3)! (a :)']: (11) 
ad 


of twist and extension are given by 





3. Twist held constant. Increasing axial load 

We have to distinguish between two different sets of critical conditions: 
one in which the twist is such that the bar is all elastic initially, the other 
in which the twist is such that the bar is plastic to a given radius initially. 
This procedure has been used for a thin tube by Hohenemser (2) to test the 
Reuss equations. 

(a) All elastic initially 

Let ¢ be the prescribed constant twist per unit length. Then, while in 
in elastic state 


oC 3G In(l h), T Grd. 


Suppose that at time to; o has increased to Oo; l to i, and that the outside 
of the cylinder is just plastic. Then 


To 3G In(I,/h), (13) 
To Gad, (14) 
md o,+375 r’. (15) 


Subsequently, when the bar is plastic to a radius c, 
oO 3G In(l/h). T Grd if < €}. (16) 


Conditions on the plastic boundary give 


347074" Y?, (17) 
For ¢ ? a, we put dé () in equation (3) and eliminate dA from (2) and 
3). Then , 
al Oo 
1G — dr+de. 
But. from (4 37 dr+oadeo 0. 
- , 6G | y-1 
Therefore on In / In| + : ! const. (18) 
) 1—oY-! 
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An element at radius r becomes plastic when / is such that 


3G In(l/h) = (Y?—3G*r?d?)! = o,, say. (19) 
Therefore from (18) 
6G { ito) fi —o. ) 
In(1// In ae 
y i(t/h) iF oY it ro, YI 
. (J \8GiF ; ‘ l oY ] o,.¥ A 
that is (7) exp(—2e,/Y) [—oF1 10, ¥"" 
Putting (1/h)6@)) B we have 
: 0, —20r/T" L (Be—2orl¥ me i 
. y] | (b 1)( Vy Y). (20) 
Y Be-2or¥ 4-1 4-(Be2or¥ —1)(¢,/Y) 


o, is given by equation (19), and then equation (20) gives ain terms of / andr. 
From (4), (18), and (19). 


: 2(G/Y rd Bier} oe 
Y Be-®r'¥ 4.1 +. (Be-2er'¥ —1)(,/¥)’ 7 

Suppose, in particular, that ¢ is such that the outside radius a is just 

about to become plastic when the twisting is completed. Thend — Y Gv3a 
and therefore 3 
Oy | 7 5 
y ( 2) 

The bar will have become completely plastic when In(//h) Then. 


at the radius a. 


ise ~ 
l—oy — 


( 


7616. 





Therefore o is within 14 per cent. of Y. If the extension is continued until 
In(l/h) Y/G. atr a 


oO 


y 


tanh 3 0-995] 


and o is within } per cent. of Y. It is therefore apparent that o rapidly 


> 


approaches its asymptotic value Y while the strains are still of an elastic 
order of magnitude. 

(b) Plastic to radius b initially 

In this case we suppose that the constant twist is such as to cause the 
bar to become plastic to radius b. We then have, prior to the application 
of the tension load, 


r Grd (7 b), T Y/v3 (r b). (22) 


Y/Gv3 ¢. (23) 


From (22) it follows that b 
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ON COMBINED TORSION AND TENSION 33 
As the axial stress is applied, 7 in the plastic region falls and the plastic 
radius moves inwards. Suppose that o is the axial stress at time ¢ and 


radius r. The differential equation is the same as that in case (a), and so 


) , TC <a 
pe In] In(; fe 4+ const. 
) l—or— 


For r a.o 0 when / h. Therefore 


F l+oY-! 


h] l—or — 
oO B ] 5) 
" Y p+) te 
On using (4) it follows that 
VOT 2B - 
aes (c<r <b). (25 
) B+] 


\n element at radius 7 becomes plastic when I is such that 


3G In (Y2— 3G?r?d2)? o 
h r 
in ease (a). Then 
Be-20r/¥ _] 1 (Be-2orlF 1 1)(¢/Y) 
r Be-2or¥ 4 14 (Be-2o¥ _1)(9/V) 
¥ : 2(G/Y )rhBte-or'¥ 
= ) Be-2or!} 1+ (Be-2o"'F¥ —1)(o, Y) 


These are exactly the same as equations (20) and (21). 
The variation of load with extension and fall of torque with extension 


an be calculated numerically if required. 


1. Fixed extension and increasing twist 

Once again we have to distinguish between two cases. In the first it is 
supposed that the cylinder is given a certain extension, which is to be kept 
constant, and which is too small to cause the bar to become plastic. The 
twist is then gradually increased from zero. In the second, the extension 


ssuch as to cause the entire cylinder to be plastic. 


(1) All elastic initially 


The axial stress is supposed to be oy initially. While the bar is completely 


elastic Y > 
Co ory T Grd. (26) 

+ The bar will become plastic to radius c during the subsequent twisting 
when 7 rT,. where 2 os ak ane 

0 V3 Tp (Y? oa), (27) 


092,17 
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Using the Reuss equations (2) and (3) we have on putting de. 0 
l do 
dx —, 
2G oO 
' * 
and Gr dd = dr. 
Y2*— 37? 


Integrating. 


G {] t NST Y 
2V3 —rob In| j+const. (car <a). 
) j j ) ( 


l—~v3r . 
When ¢ = 7,/Gr, we have 7 = 7,. and so 
),/2 ? 2 4 
-~VoO ‘ ] NOT ) ] VOT, } 
—{Grd—7,' In! ’ 0 | 
} \l—v37/Y 14+ 37,/Y} 


which gives 


V3 tanh(v3/¥Y ){Grd—r,)'+ v3 7)5/¥ (28) 
- = - \ aw 
) 1+-(v37)/Y )tanh(v3/Y {Grd—z,! 
It follows from (4), (27). and (28) that 
o To | sech(v3/¥ )(Grd To) | 
Y Y |\1+(\37,/¥)tanh(v3/Y)(Grd—7,)) 
In order to be specific, let o, = Y/2 and suppose that rd Y G: that is 
the shear strain is Y/2G. Then 
3T tanh(v3/2)+ v3/2 
: : ens Hn di 0-975. (29) 
) 1+ (vV3/2)tanh(v3/2) 


[t is evident, therefore, that 7 rapidly approaches its asymptotic value 
Y/v3, being within 3 per cent. of this value when the shear strain is only 
Y /2G. 

(b) All plastic initially 

Initially o, = Y, +r = 0 everywhere. The differential equation is the 


same as that in case (a). so that 


G ] 1 NOT r 
IV3—rdb In _+ const. 
} 1— v3 7/} 
In this case when 7 = 7). d = 0, and therefore 
. NO 
V37/} tanh — Grd (30) 
} 
Co V3 
and y sech y Grd. (31) 
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When 7¢ is only 2Y/G, 7 is 0-998Y/v3 and a is 0-0625Y. The load L and 


torque 7’ are given by 


L 
2 | Rsech RO dR, 


7Y a" 
0 








ar 1 
N37 Pee 
y 2| R*tanh RO dR, 
Tia” d 
0 
| 
0k 
} 
Q) | 
>|'e 
~ = j 
k /7 
it 1S U / 
/ 
> t i —— a SSS eee CN 
=9) 2 3 4 Sh 
Fic. 1. Variation of load and torque with angle of twist. 
alue Extension held constant. 
nly where 7'/a R and (v3G/Y)ad ®. Fig. 1 shows L/7Ya? and v3 T/7Ya5 
plotted against ®; it can be seen that they rapidly approach their asymptotic 
values zero and = respectively. 
the - . ; . 
he 5. Constant ratio of load to torque 
In this case it is supposed that o and 7 are zero initially, that the load 
ind torque are then applied, increasing from zero, in such a way that the 
ratio load to torque is constant throughout the distortion. 
Eliminating dA from the Reuss equations (2) and (3) we obtain 
1 dl lor lodr. 1 do 
(20) l dé 371 3G 7 dé 3G dé 


which, with (4). leads to 
31) 3 dl (Y2~—3r?)tr Y? ] dr 
l dé T l G +(Y2—37?) de 
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Ga 3 dl G 


We put JY =z, = 


--6, and a=—-— — 

Yh 

where « is now a varying quantity. 
The equation then becomes 


dz rh 
- = (1—3z?) — —-— a(1—32?)z. 32 
dd - l (32) 
a 3G dl 
olnce — —, 
Yl dd 
on integration we have 
3G nn") _ fae 33 
— in dd. ye 
mt (98) 
0 
While the bar is still completely elastic, 
, 2 T 
™ ‘add and = sh 
7 Ya" | as 7Ya’ i‘ 
0 
ob 
2 Xx dd 
as . ab 5 : 
Therefore — constant K, say. 
| d 
This equation is true for all ¢ in the range for which the bar is completely 
elastic; therefore « = K/2 until plasticity begins. If the outside of the bar 
becomes plastic when ¢ = dp, then, at the outside, 7/ Y = dp, and it follows 
from (33) that o/ Y vf,: on substituting in (4), we find 
xv7h2 +- 362 B (34) 


This equation determines ¢, when AK is given. 
Subsequently, when the bar is plastic to a radius c, we have 
K dp 712 9.9 
\ r vc” 9 Q~ 
; dy xdb| +— ¢* L. (35) 


¢d 


9 


The load and torque at any stage are given by 


L ¢2 ce? \h 1 
RR (1-334) "+ | (1—322)!R AR, 


7Ya* as a” ; 
T 1 /c\4 A , 
7 Ya3 ) (‘ p | 2R* dR. 
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Therefore 


1 

7, 2etla*{l—3¢c* a®)g*}!+2 [ (1-322)! RdR 

aL cla r 
T a Pee ” 
(ct/a*)d+-2 | zR?dR 


c/a 


Given K we have to determine the strain-path, that is, x. To do this a value 
of « is assumed, say for c/a = 0-9; we then determine z from the differential 























— = 4 4 4 A. 4 j 
0 0:2 0-4 06 0:8 10R 
Fic. 2. Variation of longitudinal stress as deformation 


} yrogresses. 


equation (32) and calculate a7. In general this is not equal to K and 
another « must be assumed. This process is continued until aL7'-! is 
equal to K; c/a is then taken equal to 0-8, an a assumed, and so on similarly. 
At each stage the values of z at intervals of r/a = 0-1, say, must be calcu- 
lated in the plastic zone so that the integrals in aL7'-! can be found with 
a sufficient degree of accuracy. 

In the present calculations, K is taken equal to 2, which implies « = 1 
initially and, from (34), dy 0-5. 

From section 2 we can find the limiting value of « corresponding to the 
given ratio of load to torque when the cylinder is completely plastic and 


the elastic strain-increments are negligible. This result will give some 
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where 


1 1 
r aR *" akdR 


Therefore K — : : 
re € \ (q2-+- 3 R2) } (q2-4 3 R?) 


SS 
< 
tl 

f 


a=0 
n) ee OT 
~ 10 04 02 03 0:5 % 10R 
“~!Y, 
7 a a=2 
mn . te] 
@=c 
a=] 
_0.9 
U-2 
Fic. 5. Residual longitudinal stress for different 


strain-paths. 
When K 2. this reduces to the cubic equation 
to 5a*— 12a-+4 0 
f which the relevant positive root is 0-874. 
the Figs. 2. 3. and 4 show the stresses and the strain-path respectively. 


6. Unloading, after loading with constant strain-ratios 


[It has been shown in section 2 that. for constant relative rate of extension 


and twist. the stresses in the plastic region c r <aare given by 
v) Yria 


9 


(y?-+ 3r*/a*) (a?+- 3r?/a?) 


where Y Is constant. 
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The load L and torque 7' at length / are given by equations (11) and (12), 
If the bar is now partially unloaded by removing AL and AT’, where 
0 <A <1, the residual stresses in the previously plastic region ¢ <r <a 
are given by 


o a AL R- r T R ; 2AT R 
Y (a?+ 3.R?)! Ya?’ a Y  (a?+3.R?)! Ya?’ 
0-6 | 
0:5 | a=0 
0-4 + 
0:3 | ‘unl? 
x 
y \ 
02+ \ 





| ~__a=20 --\ 
oD Of oz 03 05 Ov TOR 
\\ a=20 
YS 
—041F \ aat0 
\(aa02 
ae al a=0 
Fic. 6. Residual shear stress for different 


strain-paths. 
; o* 37" 
These stresses become negative ultimately for certain R. Should po Py 
be equal to unity, for some values of R and A, the bar will be about to become 
plastic again during the unloading. The most critical cireumstances for this 
will be when the bar was previously all plastic, that is c = 0, since then the 
residual stresses will be numerically greatest. We have, therefore, to investi- 
,; . , o 372 
gate the sign of the expression —.+ ——1 when c = 0. 
y:' y2 
done numerically for a range of values of « and for values of R between 
0 and 1, taking A 


This has been 


1. The expression is never zero for this value of A and 
can therefore never be zero for any smaller A. Hence the bar does not 
become plastic again with this simple method of unloading. 

Figs. 5 and 6 show the variation of the residual stresses with v/a for 
various values of «. 
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PROBLEMS OF CROSS-VISCOSITY 


By I. BRAUN and M. REINER 
(Rheological Laboratory, Institute of Technology, Haifa) 
[Received 9 April 1951] 


SUMMARY 

Reiner’s (1) equation of a general viscous liquid, which differs from Stokes’s 
equation by a quadratic term defining a ‘coefficient of cross-viscosity ’, is integrated 
for the cases of (i) a cylindrical rod rotating in a liquid, (ii) the rotating coaxial 
cylinder viscometer, (ili) the tube viscometer, (iv) the parallel plate-torsional visco- 
meter, thus providing criteria for crucial experiments to decide whether the ‘centri- 
petal effect in certain liquids is totally or partially due to cross-viscosity. Provided 
the liquid climbs up the rod in instrument (i), there should be radial tension upon 
the internal cylinder of instrument (ii). In instrument (iii) the wall pressure should 
be greater than the Newtonian hydrodynamic pressure by a term involving the 
ordinary shear viscosity. Finally there should be in instrument (iv) a parabolically 
distributed pressure upon the plate which vanishes at a distance from the axis of 
rotation which is 1/v3 times the radius of the plate. 


Introduction 

At the 1946 meeting of the British Rheologists’ Club, Weissenberg first 
demonstrated very strange phenomena in the hydrodynamic behaviour of 
certain very viscous liquids. These were described by Garner and Nissan 
(2) and subsequently by Weissenberg (3).+ When such a liquid is sheared 
between a rotating vessel and some coaxial inner member, which is held 
against rotation and is either fixed in position or is free to move up and 
down the axis of rotation, the liquid is drawn inwards against the action 
of the centrifugal forces and upwards against the forces of gravity, the 
whole arrangement forming a sort of ‘centripetal pump’. Weissenberg 
names such liquids ‘general’ liquids, to distinguish them from the ordinary 
viscous or ‘special’ liquid. All liquids which he demonstrated exhibited 
also elasticity of shape under the action of transient forces, and Weissen- 
berg attributes the ‘centripetal-pump effect’ to that elasticity. However, 
in a paper preceding Weissenberg’s first publication, Reiner (1) deduced 
from mathematical principles that there is a theoretical possibility that the 
most general viscous liquid might show just such behaviour. He established 
as the most general relation between the stress tensor p} and the tensor 


of rate of deformation or (better) ‘flow’ f7 defined by 


l f/ov' ev! \ 
fi 4 (0.1) 
E 2 \ C2 Cas 


Compare also Garner, Nissan, and Wood (4). 


[Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 1 (1952)] 
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PROBLEMS OF CROSS-VISCOSITY 
the rheological equation 

Pi = Kit AR+hALS. (0.2) 
where /), F,, and F, are functions of certain material constants and the 
invariants I, II, and III of f7, and he showed that the appearance of F, 
gives rise to the phenomenon of ‘cross-viscosity’, i.e. the appearance of 
stresses in directions ‘across’ the velocity gradient, or in other words to 
phenomena leading to the ‘centripetal pump’ effect. However, he con- 
cluded that this proved that F, must vanish, because such an effect ‘had 
never been observed’ (loc. cit., p.360). Actually the effect had been observed 
it about the same time by Weissenberg and others when experimenting 
with incendiary gels and was described by Russell in an unpublished thesis 
submitted to the Imperial College, of which Reiner had no knowledge. 
Later Reiner (5) showed that equation (0.2) also describes the most general 
elastic behaviour of an elastic material exhibiting an elastic effect analo- 
gous to cross-viscosity, named cross-elasticity, when ff is replaced by the 
strain ej. Many years ago Poynting (6, 7) had found that a steel wire 
lengthens on excessive torsion, a phenomenon akin to one observed later 
by Rivlin (8), viz. that a rubber cylinder which is twisted between plates 
which prevent any change of its length exerts pressures upon these plates. 
Poynting had shown that finite shear of a classical elastic solid, i.e. one 
with vanishing F, may give rise to a pressure in the direction normal to 
the shearing displacement, and this seemed to be sufficient to account for 
the lengthening of the cylinder when its ends are not restrained. One could 
therefore conclude that in an elastic liquid what was at the root of the 
centripetal pump effect was in fact the Poynting effect of its elastic com- 
ponent, but Rivlin (10) argued that elastic properties in a liquid could not 
contribute to the centripetal-pump effect and that the latter must be 


entirely due te 


viscosity as derivable from Reiner’s equation (0.2). Ina 
recent paper Oldroyd (11) has shown how from a special rheological 
equation of an elastic liquid the appearance of a centripetal-pump effect 
can in fact be derived and has added to his paper a polemical note against 
Rivlin’s view. The following problem therefore arises: derive from (0.2) 
equations describing the behaviour of a viscous liquid, devoid of all elasti- 
city of shape, in various apparatus—is it then possible to establish relations 
which may show when experimentally tested that the centripetal pump 
effect is due to cross-viscosity and not to elasticity? It is the purpose of 
the present paper to derive such relations for a number of apparatus, viz. 

(1) a cylindrical rod rotating in a liquid, 

(i) the rotating coaxial cylinder viscometer, 


Compare also Oldroyd (9). 


































44 a 





BRAUN AND M. REINER 


(ili) the tube viscometer, and finally, paran 
(iv) the parallel plate-torsional viscometer. equal 
1. Before dealing with the theory of these apparatus it is instructive to | the 
treat the case of simple laminar shear between a stationary platen (y = 0) | bY ™ 
and another parallel to the first at a distance D, moving in the directiony, 4") 


with the velocity V, so that the velocity gradient ment 
; vlin 
G = dvu/dy = V/D. (1.1) | * Ty} 
When a steady state of flow is reached, the velocity is elent 
v, = Gy, v, = v, = 0, (1.2) — displ 
and the flow tensor is, in accordance with (0.1), In di 
, mixe 
0 1 0 

G | elast 

fi=s\}1 0 Of (1.3) 
“lo 00 || cale' 
, are | 
Introducing the expression (1.3) into the equation (0.2) we get _ 
gl® 1 Ol gilt 9 off com 
pi = KH8+F>|1 0 OF 79 1 Of. (1.4) whe 
~|0 0 0 }0 0 0| evli 


We now have to introduce these expressions into the stress equations | lle 
and integrate. When doing so we must keep in mind that while F, and F, | hall 


can be considered as material ‘constants’, this is not the case with regard 2 
to F, because this must reduce for vanishing rate of deformation to the vist 
‘hydrostatic’ pressure.} /, must therefore generally be a function of the lon, 
coordinates. In our case we find from the third stress-equation, neglecting uni 
body-forces, that p,, is not a function of z. If the liquid is exposed to air | the 
pressure —II at some end z = z,, the boundary conditions then yield mo 
Pe, = F, i. (1.5) | inf 

- ‘ , : . se 
Considering this, (1.4) now gives i 
We 


Dae Pus F, G?/4—T] | 


7 y ) i (1.6) on 
Pry Pyx F, G/2 


Tm 
Let 7 be the force which must be exerted upon the platens of area A in 


the direction x in order to maintain the movement, then evidently pu 

T = AF,G/2, (1.7) | 38 
and therefore F, = 2n, (1.8) cs 
where 7 is the ordinary shear-viscosity of the liquid as in the ‘special’ to 


Newton-Stokes viscous liquid. In addition, however, the appearance of the fu 
— : a in 
+ For the difference between hydrostatic or thermodynamic pressure and mean pressure 

Pm = Paq/3 compare Truesdell (12). ) th 





Live to 


tion a 


(1.2 


«) 


(1.4) 


tions 
nd F, 
card 
0 the 
f the 
cting 
O air 


(1.5) 


(1.6) 


A in 


ssure 





PROBLEMS OF CROSS-VISCOSITY 45 


parameter F, carries with it stresses ‘across’ the direction of shear which are 
equal in both (x and y) directions of the cross-section. This is generally not 
the case with cross-elasticity (compare Reiner (5)). One could measure p,,, 
by inserting suitable pressure gauges into the platens, but it is difficult to 
carry out actual experiments on these lines. The experimental arrangements 
mentioned above have therefore been devised in viscometry involving 
cylindrical and polar coordinates. 

The expressions in these coordinates are well known in the case of 
elasticity; compare, for example, Love (13) or Reiner (14). Replacing the 
displacement ‘wu’ in these expressions by ‘v’, we can use them in our case. 
In doing so we must keep one point in mind. Equation (0.2) is written in 
mixed tensor-components, while the equations of the classical theory of 
elasticity or hydrodynamics use so-called ‘physical’ components. In the 
calculations of the present section, where Cartesian rectangular coordinates 
are used, this difference is immaterial. It has, however, been pointed out 
by one of us (Braun, 15) that we can replace in any equation the mixed 
components of a tensor by its physical components, and vice versa also, 
when curvilinear but orthogonal coordinates are used. This is the case with 
cylindrical and polar coordinates which we are going to use. Following 
Ollendorf (16), we shall denote physical components by placing indices 
half-way up and putting them in brackets as, for example, pcr). 

2. We consider first the problem which Newton treated for a simple 
viscous liquid in his Principia, Lib. ii, Sect. ix, namely ‘If a solid infinitely 
long cylinder in a uniform infinite liquid revolve round its axis with 
uniform motion, and the liquid be caused to revolve by the impulse 
thereof only, but every part of the liquid persevere uniformly in its 
motion....’. Newton neglected mass forces and we shall do likewise. The 
influence of mass forces, gravitation, and centrifugal forces will be treated 
separately, but it is immaterial for the purpose of our present investigation. 
We shall, however, in contradistinction to Newton, take into account not 
only shear viscosity but also cross-viscosity and our liquid is accordingly 
non-Newtonian’. 

We are considering this problem in the first place because the centripetal- 
pump effect appears here in its most striking form. When a cylindrical rod 
is rotated in the liquid showing such an effect, the liquid will climb up the 
rod. To simplify calculations, let us assume that the liquid is contained 
in a coaxial cylindrical vessel so wide that it can be assumed as extending 
to ‘infinity’. The rod will, of course, not be ‘infinitely long’. Let us 
furthermore assume that there is a lid at the surface of the liquid cover- 
ing it except at a narrow annular slit at the edge of the vessel. All 
these assumptions will introduce ‘end effects’ which for the sake of 
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simplicity we neglect, but which can be dealt with by the usual 
methods. 
Applying a ‘semi-inverse’ method, we assume 
v, 0. v9 Tw, v, = 0, (23 


where w is a function of r and z only. 


0 1 0 
The flow tensor is fi Allil 0 O (2.9) 
lo 0 0 
Cw 
where A steal (2.3) 
” OF 
1 0 0O 
We have po iy Asii0 1 Of. (2.4) 
0 0 0 
Equation (0.2) now yields 
f+, A? 2A 0 | 
DP} 2n | F | f A2 QO | (2.5) 
0 0 F, 
The stress equations ((11) in Reiner (14), p. 41, or Love (13), p. 91) become 
oF,/ér + F, eA?/er ° | 
C F, oo-+ 2nr cA /er- 4A i) 


: (2.6) 
ok, Cz 0 | 


Considerations of symmetry require that @F,/e0 = 0 and the second 
equation accordingly vields 
w C;"/r?+C,. (2.7) 


The boundary conditions are w = 0 for r = «©, w 


I~ 
Ww 

a 
ow 

e+ 


radius of the rod. This makes 
w = OR?*/r? (2.8) 
as in a Newtonian liquid (cf. Reiner (14), p. 26). The first and third of 
equations (2.6) now vield 
F, = Pu F, OQ? R4/r4+ K. (2.9) 
This confirms that F, will generally be a function of the coordinates. 
Because of the slit in the lid the liquid is there in contact with air of which 
the pressure is —II]. Therefore p,. Il for r = «© and 
[1 +-p,. F,Q2R4/r4, (2.10) 
If F, is assumed to be positive, the lid exerts a pressure upon the liquid 
which at every point is given by (2.10), vanishing at infinity and largest 


near the rod. Therefore, if we remove the lid. the liquid must rise near the 
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rod. We accordingly get the climbing-up effect. It would be difficult and 


usual 
the effort not warranted to try to calculate the free surface. When the liquid 
rises, gravity will at the same time cause it to flow down. The combination 
(2.1) | of both tendencies results in a circulatory movement, the liquid rising at 
r the side of the rod and flowing down at the free surface away from the rod. 
[his is complicated by the effect of centrifugal forces which cause a 
lepression around the rod. The main point is this: if the climbing up is 
2.2) due to the appearance of F,, the latter must have the positive sign. This 
fixes the sign of F, for conditions in other apparatus. 
(2.2 3. Proceeding to the rotating coaxial cylinder-apparatus of the Couette- 
. Hatschek type, the calculations for the rotating rod which showed us that 
the kinematics of the case are the same as for a Newtonian liquid suggest 
(2.4 that similar conditions may exist here. Let R; be the radius of the internal 
stationary cylinder, #, the radius of the external cylinder rotating with 
elocity Q, then the only non-vanishing velocity component is vg = rw, 
where w = Q(1/ R?—1/r?)/(1/R?—1/R?), (3.1) 
2.5 in accordance with Reiner (14), p. 28. 
We need not fear to have made a wrong assumption. If we calculate the 
‘ome | stress tensor on the basis of (3.1), we shall be checked by the requirement 
f compliance with the stress-equations. 
(26 The flow tensor is as in (2.2), with 
A Q)/r?(1/ R?—1/ R2), (3.2) 
ond f+ 5, A? 2nA 0 
from which pi 27nA F\+F, A? 0 |. (3.3) 
2.7 . : a 
the With this, the stress equations are 
oF, /er+ F, dA2/dr = 0) 
(2.8 oF. a0 = —2ylr*4 4 2A) 0 | (3.4) 
d of dr 
C F, Cz 0 
2.9) [he second equation confirms (3.1), the other two yield 
> 2 2 2 
se p Fy x ee a+ K. (3.5) 
hich iid | "So R?)2 - 
Here again we assume the liquid covered by a lid attached to the stationary 
10) ylinder and leaving a slit open to the air near r R,. We shall find 
juid Il for ? R, and therefore 
est R40O2R 


K . en ame (3.6) 
the | (Rt Rey 


i 
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RIQ2F, 
(R2— R?)2 


As before, the lid therefore presses down the liquid near the internal 


This makes Il +-p,. (R4/r4#—1). (3.7) 


cylinder, and if we remove the lid the liquid will tend to climb up on the 
internal cylinder. Such an effect was observed years ago by Dr. Lee (oral 
communication) when testing certain bitumens in his coni-cylindrical 
viscometer (Lee and Warren (17)), an effect which was considered as a 
‘nuisance’ in viscosity measurements. 

We now calculate por) = F)+ F, A? and find 

QQ? Rt 
(R2— R?)? 


This is a tension, but we cannot be sure that its magnitude is correct because 


I] + porr) F.. (3.8) 


if we assume a lid with slit near the internal cylinder, or p,. = 0 forr = R,, 
we find Q2R4 


I+ pr) F.. (3.8.1) 


(R2— R?)? * 
This is also a tension, but larger in the proportion of R4/R*. When there 
is no lid the radial tension will have a magnitude 
Q)? 
(R?2— R?)? 


’ 


fi-j p(rr) F,aR}, (3.8.2) 


where a@ varies between | and (R,/R;)*. 

If one could measure pcrr) for different radii R; and R,, one could deter 
mine a and F,, and at the same time check the theory. A modified Couette 
apparatus has been designed by Dr. Frei and built at the Weitzman 
Institute of Science at Rehovoth, details of which will be published else- 
where. Its essential feature consists in openings in the inner cylinder which 
is filled with the liquid under test. At the same time the inner cylinder is 
covered by a rigid lid in which a pressure gauge is inserted. The gauge 
accordingly measures pirr). Therefore, should the centripetal pump effect 
be due to cross-viscosity, the said gauge should show a vacuum. This may 
provide a crucial experiment to decide between cross-viscosity and cross- 
elasticity. 

4. We now consider conditions in the tube viscometer. To exclude the 
action of gravity, let the tube of radius R and length / be horizontal, 
connected with a reservoir in which a pressure head —p, (including air 


pressure) is maintained and emptying to air pressure, whichis taken as —II. 


Laminar flow will be telescopic with velocity v = v(r), independent of z 
‘ ). The ris 
ind @. The flow i: 001 
fi=All9o 0 Oj, (4.1) 
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where A y’ /2, (4.2) 
l 0 0 

from which pat A*||0 O O}f. (4.3) 
0 ® .3 


Equation (0.2) now yields 


Pi 0 FR oO ij. (4.4) 
2nA 0 F- 
[he second of the stress equations then shows that F, is independent of @. 


lhe other two yield 


v C, Inr+C,r?/2+-C,. (4.5) 

As » must be finite for 7 0, C, must vanish. Assuming that the liquid 
idheres at the wall (r R), we get 

v C( R*—r?)/2. (4.6) 


The velocity distribution is parabolic as in the case of a Newtonian 
iquid and is therefore not affected by cross-viscosity. 
Introducing v from (4.6) into the stress-equations and integrating vields 
F) 2Cnz—3Cr?/8.F,+K. (4.7) 
With this, p., becomes 
Pp. 2C'nz—C?r?/8.FL+K. (4.8) 


The constants C' and A can be determined from 


| p..2rz di R?rp, for z 0 and 11 R22 for z l, 


f 0 
ielding " » 
( (Po {| 2yl | (4.9) 
K Do+ F, R2(py—I1)?/ 647222 J 
nd therefore 
F, p (l—z)/l+- ( R?—6r?)/64721?. F(p y—1)?—Mez/l. (4.10) 


Let Qt be the quantity of liquid which passes the cross-section of the 


tube in unit time, then from 


v p I1)( R?—r?) 4y/ (4.11) 


0 


we find Q/t = (py—T1) RI /8l, (4.12) 


nd we see that Poiseuille’s equation is not affected. In order to determine 


the coefficient of cross-viscosity, we insert into the wall of the tube 


ressure-gauges and thereby measure p,, at r R. But 
Pp Pp (1 m) l Rk? 2r*) 6477/7. Fy( po 11)? Iz / (4.13) 
nd prr),7 py(l—z R*( pp —l1)?/64n71?. Fy—M2/1. (4.14) 


E 
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The first term on the right side is the pressure gradient of a Newtonian 
liquid and we see that to this is added — R?(p y—II)?/647/*. FF. As F, is 
positive, the second term is a pressure, as is the first term. 

When the liquid emerges from the tube, i.e. for z = /, the first term 
vanishes but not the second. The result must be a ‘swelling’ of the 
emergent column of liquid, as has been observed by Merrington (18) in 
the case of rubber solutions. It is true that Merrington attributes the 
phenomenon to the elasticity of the solution, but we have shown that such 
an effect might be due to the cross-viscosity of an otherwise inelastic 
liquid. We meet here the same situation as in the other cases treated: 
as long as the observations are qualitative only they might be equally 
attributed to elasticity or cross-viscosity. Only quantitative measurements 
can decide between the two possibilities. Here conditions in the tube 
viscometer are particularly interesting. Not only would it be possible to 
vary R and/or /, but by determining the temperature-dependence of F, 
and of » in the coaxial rotating cylinder viscometer one could check the 
‘cross-viscosity pressure’ in (4.14), which thereby provides another crucial 
experiment. 

5. Finally we consider the parallel plate viscometer in which the liquid 
is contained between a plate rotating around a vertical axis with angular 
velocity Q and another parallel one at a distance H which is stationary. 


We assume 5 
Vg rw, (5.1) 


where w is a function of z only, while v, and v, vanish. Then 


0 0 0 
f; rem oO 21, (5.2) 
oe t.9 
rdw = 
where A : (5.3) 
2 dz 
0 0 O 
from which Rii= Ate 1 Ot. (5.4) 
a 
Ep 0 0 
We now have pi 0 F+A2% 2nA 7 (5.5) 


0) 27nA K+ A2F, 
From the second stress-equation we find 
Ww Ol —2/ff). (5.6) 


The other two equations yield 


pir = Fy FE 
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At the edge of the plates where r = R, the liquid has a free surface, more 


or less vertical, where pirr) can be put as equal to —II. This makes 


. R222 _ aoe 
K WE F,—tl (5.8) 
.] Q)? > 9° ’ ~ 
an I{ + FP . y/o r)F, (5.9) 
02 : ner 
therefore Il +-p.. F,+A?2F, spp 3r?)F,. (5.10) 


This shows that p., is distributed parabolically over the width of the plate 


with pressure at the centre, changing to tension atr = R/v3. An apparatus 
was built at this institute with a number of standpipes inserted in the upper 
stationary plate and the pressures thus determined. The zero-pressure 
point was clearly noticeable. The experimental observations will be pub- 


lished elsewhere. 


6. We may now reply to the question put at the end of the introduction 
is follows: 

Should the peculiar properties of Weissenberg’s ‘general liquids’ which 
we comprised in the word ‘centripetal-pump effect’ be due to a second-order 
viscosity term, quantitative observations in three different kinds of visco- 
meters can provide crucial experiments. These are: 

a) The torsional parallel plate instrument which permits measurement 
of the pressure exerted upon the stationary plate. The coefficient of cross- 

iscosity F, can be determined from the pressure at the centre p 


SH* 


2 22¢)2?max> 


max by 


leans of 


where H is the distance between the plates, R their radius, and Q the angular 
elocity of rotation. If the latter is not too great, so that centrifugal forces 
in be neglected, the pressure should be parabolically distributed over the 
diameter of the plate, vanishing at r R/\3. By varying the temperature 
it which the experiment is performed, F, can be expressed as a function of 
the temperature in an empirical manner (in the form of a table or a graph). 
b) The value of F, thus calculated can be checked against measurements 
radial stress exerted upon the stationary internal cylinder of a rotating 
-axial cylinder instrument. This should be a tension of magnitude (apart 


irom atmospheric pressure) 





‘e ais a factor 
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(c) Ina tube viscometer which permits measurements of wall pressures, 
these should be greater than the Newtonian hydrodynamic pressures (apart 
from atmospheric pressure) by R?p§/6477/".F,, where py is the driving 
pressure and 7 the viscosity. If the dependence of 7 on the temperature 
is known, such measurements, varying temperature as a parameter, should 
provide the final crucial experiment. It should be noted that we have 
assumed that F, and 7 are constants; the possibility that both F, and 
are functions of the second invariant of {7 is, however, not excluded. Both 
will naturally also depend upon IT. 

From preliminary observations which we have made with rubber-toluene 
solutions, it seems that the phenomena under discussion have both an 
elastic and a viscous component. It seems that the elastic component is 
more pronounced but less stable, gradually disappearing in continuous 
shear and restored through rest: a kind of elastic thixotropy. The viscous 
component is apparently stable, but much less pronounced so that it can 
be overlooked, as happened to us (compare Braun, Frei, and Reiner (19)). 
In the rotating cylinder apparatus the elastic components seem to cause 
pressure upon the internal cylinder. After some time, due to its destruction 
through shear, the pressure disappears and changes into tension (oral 
communication by Dr. Frei), which can be mistaken for the result of 
centrifugal forces but actually seems to be an indication of cross-viscosity. 
It therefore seems to us that a clear understanding of the phenomenology 
of these effects will not be gained before the manifestations of cross- 
elasticity in the above-mentioned apparatus are analysed analogously to 
those of cross-viscosity as presented in the present paper. We have such 
an investigation in hand. Only after both elastic and viscous components 
are separately and quantitatively comprehended will it be possible to 
advance a theory of the structure of the liquids to which the phenomena 
might be due. R. 8. Rivlin (20) has proposed a theory of solutions of 
macro-molecules in which he attributes the coefficient of cross-viscosity to 
the orientation of the molecules in the flowing solution. Experiments on 
the lines proposed in the present paper would provide a means of checking 
this theory. 
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NOTE 
We are grateful to Mr. Carmi for drawing our attention to the fact that —p,—II is not 
tal pressure required to maintain a continuous flow of liquid in the tube while emptying 
reservoir. An additional pressure head is of course required for imparting to the liquid, 
vhich in the reservoir is at rest, the kinetic energy of flow; this forms the well-known kinetic 
nergy correction. In addition, in our case, the liquid is compressed when forced to enter 
tube and this requires an additional head. The total reservoir head will therefore exceed 


radial pressure at the wall of the tube 
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SUMMARY 

In this paper a theory of flow at high Mach number past thin aerofoils is developed 
and a critical survey of previous papers given. The essential feature of the paper is 
the establishment of a general similitude, given first by Hayes, which expresses the 
fact that the problem of determining the two-dimensional flow at high Mach number 
past a thin aerofoil is in the limit analytically identical with that of unsteady one- 
dimensional flow produced by the motion of a piston. The error involved here is 
shown to be a factor of the form 1+ O0(M,?), where ™M, is the undisturbed stream 
Mach number. Using this similitude a note on the effect of entropy is given. 


1. Introduction 
THE present paper deals with two-dimensional rotational gas flow at high 
Mach number past thin aerofoils, with the proviso that the Mach numbers 
are not so large as to cause the melting of any bodies placed in the flow. 
For this type of flow we shall refrain, unlike some authors, from using the 
term ‘hypersonic’, because a discussion on the behaviour of rarefied gas 
flow, where the mean free path is comparable with the aerofoil chord, is 
not intended, although the condition that this effect should be negligible 
is investigated. Our attentions are therefore confined to gas flows at high 
Mach number for which the undisturbed stream Mach number, M,, assumes 
values of about 5. Besides this numerical estimate of the Mach numbers 
considered, it must also be remembered that the maximum stream 
deflexion, 5,,, of the wing is limited, being bounded above by considerations 
of drag and below by considerations of strength, so that values of about 
5-10° are taken as typical of the types of aerofoils employed for the 
application of the theory to be developed. It is seen from the above 
numerical considerations that, for this type of flow at high Mach number, 
the product My 4,., where 6,, is now expressed in radians, may be reasonably 
assumed to be of the order of 1. 

The study of hypersonic flow has attracted the attention of several 
authors, earlier ones being content to study the extreme case of an infinite 
value of the undisturbed stream Mach number. von Karman (1), by 


considering y, the ratio of the specific heats, to take the value 1, has 


[Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 1 (1952)1 
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pointed to an analogy between the dynamics of hypersonic flow and 
H Newton’s corpuscular theory of aerodynamies,t finding an approximation 
for the difference between the pressure at any point on the wing surface 
ind the undisturbed stream pressure as a constant multiple of the squart 
of the angle of inclination, instead of the linear law of Ackeret. Sanger (2) 
has used this concept to design the best wing and body shapes for flight 
it infinite Mach number. Such a method of approach, however, is not 
satisfactory in forming the basis of a reasonable approximate treatment for 
flows where the Mach number ™, is large, but M,4,, is not. 
Another approach to the problem has been given by Tsien (3), who 
deduces the similarity laws of hypersonic flow by a direct study of the 





= 
te ‘ equations of motion for irrotational gas flow. He uses the fact that for 
s the hypersonic flow over a wing the disturbed flow due to the presence of the 
mber | wing is confined within a narrow region close to the body, ‘the hypersonic 
o" } boundary layer’. The thickness ratio, 8,, of this layer can be shown to be 
oar of the order of 1/M,. An expansion of the coordinate normal to the 
| direction of the undisturbed stream is therefore valid, as it is in the similar 
| case of ordinary viscous boundary layer theory. Applying this expansion, 
the equations of motion for irrotational flow are reduced to a single non- 
high dimensional equation depending on a single parameter M,4,,. Tsien is then 
ber led to state the following similarity law: If a series of bodies having the 
low same thickness distribution but different maximum stream deflexions, 6,, 
‘the | are placed in different flows such that the parameter M,6,, remains con- 
gas | stant, then the flow patterns are similar. With this law as basis, Linnell (4) 
d, is investigates the flow at high Mach number through a straight oblique shock 
rible ind the Prandtl-Meyer expansive flow around a smooth surface and 
high deduces expressions for the force and moment coefficients for a sharp 
mes edged two-dimensional aerofoil with particular reference to a flat plate, 
bers single and double wedges, and a parabolic profile. The flow, behind the 
eam shock occurring at the leading edge, over the wing surface is, however, 
ions issumed isentropic. This assumption, particularly in the case of the 
out parabolic profile, may not be sufficiently accurate in the case of flow at 
the irge Mach number, where strong shocks are contemplated (the shock 
ove strength is shown to be of the order of M,4,,). Lighthill (5) has shown that 
ber, for two-dimensional hypersonic rotational gas flow produced by the passage 
bly fa uniform hypersonic stream through a stationary shock of nearly 
iniform strength, the reflection of the disturbance waves due to any 
eral urvature of the aerofoil surface at the shock plays an important role. 
nite However, by taking the boundary of the isentropic flow to be at the shock, 
by { Linnell deduces expressions for the pressure over the aerofoil surface, 
has t See Principia, lib. ii. sect. vii, prop. 34. 
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which are true neglecting the fourth power of the shock strength and which 
take into account the influence of the known maximum entropy change, Si), 
across the shock, which occurs at the nose of the aerofoil. The influence 
of this entropy change is not studied by Linnell, and the present paper, 
although not seeking a higher approximation than Linnell’s, lends quanti- 
tative support by the numerical study of an alternative expression (6) for 
the pressure over the aerofoil surface which shows in a clearer manner the 
effect of entropy. The ratio of the pressure at any point of the wing surface, 
given by the simple wave theory, to the undisturbed stream pressure is 
modified by a factor of the form 


where c,, is the specific heat of the gas at constant volume. For y = 1-4 
and M,6,, = 0-8 this expression becomes {1+ 0-01}. The effect of the 
known maximum entropy change is therefore small; hence the expressions 
derived by Linnell form a reasonable approximation for the values of the 
parameter MM, 6... 

A very interesting note on the subject of flow at high Mach number was 
struck by Hayes (7), when he expressed the fact that the problem of 
determining the two-dimensional flow at high Mach number past a thin 
aerofoil is, in the limit (when MM), but not M)4,,, is large), analytically 
identical with that of unsteady one-dimensional flow. Hayes extends this 
concept to flow at high Mach number past any slender body and also 
mentions that this result holds even when shock waves and the resultant 
entropy changes are present. Shen (8) has used this more general similitude 
to carry out numerical computations dealing with the flow at high Mach 
number over a slender cone. The present paper, by a more detailed study 
of the exact equations of motion, establishes this approximate similitude 
for two-dimensional flow at high Mach number past a thin aerofoil and 
examines the errors involved. The actual error involved is found to be a 
factor of the form 14+-O(./,*); hence such an approximation will be of 
great value. From this similitude the expressions derived by Linnell are 
obtained from the unsteady one-dimensional theory, considering the 
variations of entropy from the known maximum entropy change at the 
shock to be small. 

Before entering into a discussion of the equations of motion it is desirable 
to study first the physics of the problem, in relation to the magnitude of 
the mean free path and the effect of viscosity. If the mean free path is 


small in comparison with the dimension of the flow field, which in our case 


is taken to be the thickness of the viscous boundary layer, then the gas 
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vi 
can be considered as a continuous medium and the theory of conventional 
cas dynamics applied. The case for which this is not so will not be treated 
here, since considerable modification is required and the field of super- 
aerodynamics must be entered (9). Also, since the region of disturbed flow 
is thin, the presence of the viscous boundary layer, and hence of the 
viscous terms in the equations of motion, may be felt. The resulting 
influence of these two factors is somewhat linked together, for it is now 
shown that if one can be neglected so can the other. 

For a body of length L placed in a stream, which has velocity U and 
kinematic viscosity vy (= py/ Py Where po is the coefficient of viscosity and 
po, is the density), the appropriate Reynolds number, R&,, is defined as U L/v. 
For large Reynolds number it is well known that the thickness ratio of the 


viscous boundary layer, 


5, = 1/VR,. (1.1) 
Thus, since the thickness ratio of the region of disturbed flow (the hyper- 
sonic boundary layer) is approximately equal to 1/./,, we may write 
©... M 
=. (1.2) 
0, VR, 
Now, considering the molecules to be rigid perfectly elastic spheres, the 
mean free path can be calculated (10) as 
l L-255,/y(v.M, U), (1.3) 
which, when combined with equation (1.1), gives 


l My 


~ 


~ ; 
LS. VR 


t ( 


(1.4) 


Comparison of equations (1.2) and (1.4) shows that if R, M?, then 
neither the presence of the viscous boundary layer nor the magnitude of 
the mean free path has an appreciable effect on the dynamics of the flow. 
In particular, only at extremely high altitudes will the mean free path be 
comparable to the dimension of the flow field. At an altitude of 50 miles 
Maris (11) calculates the mean free path to be nearly 1 in., and hence the 
lir there certainly cannot be regarded as a continuous medium. 

The author wishes to express his thanks to Professor M. J. Lighthill for 


his continued help and encouragement in the preparation of this paper. 


2. Boundary conditions 

In order to ascertain approximate orders of magnitude for the flow 
quantities in the region of disturbed flow, it is necessary first to study their 
magnitude at the surface of the aerofoil and at the shock. It is then 
reasonable to assume that these orders of magnitude hold throughout the 


region of disturbed flow. 
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Consider an aerofoil, whose profile is given by y = y,,(2) with reference 
to a suitably chosen Cartesian coordinate system such that the x-axis is 
along the chord of the profile, placed at zero incidence to a uniform stream 
which has velocity (U,0), density po, pressure pp», and appropriate Mach 
number M,. Let the velocity, density, and pressure in the region of 
disturbed flow be (U +-u, v), p, and p, respectively, where u and v are small 
in comparison with 

The condition for zero normal velocity at the upper surface of the 


aerofoil is 
dy,,, 


ey 
alee 


0. (2.1) 
Now since w is small and y,, is of the order of 5,,, it is seen that v = O(U8,,). 
This order of magnitude or smaller is assumed to hold throughout the 
region of disturbed flow and in particular at the shock. 

At the shock, y = y,(x), the Rankine shock equations can be written as 


dy. dy. 
(py— p)U Ss of / ) (2.2) 
aa Aa 
dy. 
0 i, 2.3 
Uu a i ( ) 


as dy,\? ,)f dy,\*) 
U2 1 — Po) (H D2 24 
po U2(1—P) (TR) = (n—pol(1+ (7) | ea 
Uu- ku? | 932) Y f P| (2.5) 
Y 


By process of elimination equations (2.2)-(2.5) reduce to a_ single 
quadratic equation, from which the slope of the shock can be given 
explicitly in terms of the angle through which the flow is deflected on 
passing through the shock and initial Mach number, the root appropriate 
for flow at high Mach number being 


dy, Yr lw a | yA 1G ;) + 1 | (2 6) 
dx 4U'A\\ 4, (7; M2} a 


Hence dy,/dx is of the order ‘ l | Mo, and noting, from the definition of the 
Mach number, that py = O(p,) U?/M?), we are led from the above equations 


to the following table of orders of magnitude 


TABLE 1 


Quantity v Ye u Pp 


Order of magnitude US, 1/M, Us 
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Using non-dimensional flow quantities given by 


v U5... 0’, u (U6,,./.M,)u’, Po (py U?/ M2) po, 
Pp (Po l "2 M?)p’, (2.7) 
together with x Dx’, y, = (L/M5)y;, (2.8) 
equations (2.2)—(2.5) become 
(1 Po l dy. - 0 # (2.9) 
p M, 4,, dx M* 
dy’, 
0 u’ +-9' =. 2.10 
dx’ ( 


) dy’ : ’ l 
es) j p’— po o( ) (2.11) 
| p (—" Ms 


w’ 1»; 2 Y I __ Po , ! O l 5) 12 
M3, *” ~ y—1 URS (x Pi) _- ee 


Hence, neglecting terms of O’1//?) we find that 


2% ts = : 
Y * iz u) 1] /: M5, ", (2.13) 
)’ 2» [dy’\2 vy l 
P rol . (2.14) 


It is seen from the above equations that the similarity laws still hold at 
he shock, the flow quantities depending on a parameter M, 6,,. 


In terms of the original flow quantities, expressions (2.13) and (2.14) be- 


come 
v 2 | di d 
( J" . ; Ys (2.15) 
y+1}\dx MP dx 
p “) wee) Y L (2.16) 
Po ; | dx y+ ] 


Expressions (2.15) and (2.16) correspond exactly to those obtained by 
Linnell (4), but the method of derivation given here is superior in the sense 
that the errors involved are known. From equations (2.15) and (2.16) and 
lable 1 it is seen that shock strength ((p—po)/po) is of the order of M, 5, 

\ssuming that the order of magnitude for uw given in Table 1 holds at 
the upper surface of the aerofoil, the boundary condition expressed by 


M2), 


equation (2.1) now becomes, on neglecting terms of O(1 
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3. The equations of motion 
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Since we contemplate the formation of shocks whose strengths may not 
be so particularly small for the neglect of entropy changes to be regarded 
as a reasonably good approximation, the equations of motion for two- 


dimensional steady inviscid rotational gas flow are considered: 


(wwe? ca ou mh ‘ (the equation of (3.1) 
Cx cy Ox oy continuity); 
4404-2 =o, (3.2) 
On cy pox (the equations of 
; , ie momentum); 
iw 140-4. nt (3.3) 
Cx Cy poy 
together with the equation of state for a perfect gas 
F — PO eSie, (3.4) 
py Po 


and the equation expressing the assumption that the changes of state are 
adiabatic 
—+tWv v: (3.5) 
Ox oy 
here S is the entropy increase in the disturbed state and c,, is the specific 
heat of the gas at constant volume. 

It is now assumed that the approximate orders of magnitude given in 
Table 1 hold throughout the region of disturbed flow, so that substitutions 


in terms of non-dimensional quantities given by equations (2.7) and (2.8) 


remain valid. This leads to the following non-dimensional equations of 


motion: _* ‘ ay’ 
eed aa A a (3.6) 

M,5,, 0a oy’ — * Oy’ | Ms 
a rd . +e - ~ +O 4 = 0, (3.7) 

M5, ox’ © oy’ p M282. ox M* 
L ” e . Po a = LO J _ 0. (3.8) 

My 5,,. x cy’ pp Mid. oy M? 

oS S S 
5 te’ — 4-Gl | oe &, (3.9) 

M,4,,, ex oy M; 


Neglecting terms of O(1/.M?) and assuming that M,6,, = O(1), we ob- 
serve that the above equations are characterized by the parameter M)4,,, 
thus indicating that the similarity law for flow at high Mach number 
developed by Tsien is still true, even when the effect of entropy is taken 
into account. 
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tewriting in terms of original flow quantities and neglecting terms of 


> 


O(1/.M?). equations (3.6), (3.8), (3.9) become 


UF 4 2 i 0, (3.10) 
Ca ct C 
( cl l op ‘ 
| v 0 (3.11) 
¢ cy p cy 
cs os 
1 a 0 (3.12) 
« cy 


which, on writing a ('t, correspond exactly to the equations of motion 
governing one-dimensional unsteady anisentropic inviscid gas flow in the 
direction of the y-axis. By the same transformation, the boundary con- 
dition (2.17) at the upper surface of the aerofoil corresponds to that 
holding at a piston moving with velocity proportional to the slope of the 
verofoil: and the boundary conditions (2.15) and (2.16) correspond to those 
ippertaining in one-dimensional theory, at a normal shock moving with 
. velocity which is equal to Ll’ times the slope of the shock in the two- 
dimensional theory. 

Hence. neglecting terms of O(1 


Me), 


holds for anisentropic flow in two dimensions and can be expressed thus: 


Hayes’s more general similitude 


If an aerofoil, whose profile is given by y = y,,(x), is placed in gas flow 


it high Mach number, then the problem is analytically identical with the 


unsteady one-dimensional flow problem in y and ¢, where 


f ai (3.13) 
resulting from the motion of a piston with velocity U dy,./dx at time t, and 
the problem is characterized by the parameter ™, 5... 


Thus the problem of two-dimensional flow at very high Mach number 


c 


past a thin aerofoil has resolved itself into the study of unsteady one- 


flow. about which there is a considerable literature. The 
of the form 1-4 O(M, 2). 


here that, once v and the velocity of sound a 


dimensional ga 
error involved here is a facto 
It is worth while to mention 


, (yp p); have been determined, wu may be immediately calculated from 


the approximate energy equation, which can be determined on neglecting 


terms of O(1/.W5), from equations (3.1)—(3.5) as 
lx : (a2 —a?), (3.14) 
y—1 
4. A note on the effect of entropy 
Mention has already been made of the fact that previous authors writing 


on the subject of flow at high Mach number have given little quantitative 


support to the assumption made that the flow behind the shock is isentropic, 
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so that the truth of its validity remained open to discussion, since in gas 
flows where the Mach numbers are very high strong shocks are contem- 
plated. 

First, let us consider the numerical values of the maximum change in 
flow quantities (these are denoted by suffix 1) on passing through the 
shock. For an aerofoil causing a maximum stream deflexion of 10° in 
air flows having different undisturbed stream Mach numbers, the results 
of section 2 and equation (3.4) lead to the following table of approximate 
values for the changes in flow quantities on passing through the shock. 


TABLE 2 


0 1 (Py o)/ Po O/C (Sy Sq)/l (y¥1—Y0)/l 
c o'16 1°75 0°043 0°0072 O*1O072 
o*160 2°33 0064 0°0094 o'1694 
7 0°16 5 0-098 0°01 27 0°1727 


Columns 5 and 6 represent the changes in the Riemann invariants (refer- 
ring to unsteady one-dimensional theory) defined by 


y i+: 8 —-. (4.1) 


> 
= 


In working out the pressure over the aerofoil surface, Linnell takes into 
account the maximum change in the entropy, S,, and the change in the 
Riemann invariant s,; the variations from these changes are neglected. 
The above table, however, shows that the maximum changes, themselves, 
are quite small and therefore we may expect the variations from these 
maximum changes to be insignificant near the aerofoil surface. In fact, 
the effect of the maximum entropy change on the pressure over the aero- 
foil surface is to modify the ratio of this pressure, given by the simple wave 


theory, to the undisturbed stream pressure by a factor (6) of the form 


{ 4 (3y 5)S, ; (4.2) 
8(y—1)ec,) 
which for y 1-4 and M,4,. = 0-8 becomes {1+-0-01}. Hence the effect 


of the maximum entropy change is already small, and therefore, the 
expressions, derived by Linnell, for the pressure over the aerofoil surface, 
form a reasonable approximation. 
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THE UNSTEADY, SUBSONIC MOTION OF A SPHERE 
IN A COMPRESSIBLE INVISCID FLUID 
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SUMMARY 

In this paper we investigate the effect of compressibility on the work required to 
accelerate a rigid sphere, immersed in a fluid, from rest to a final uniform velocity. 
Viscosity is neglected, but it is shown that for the accelerations for which com- 
pressibility is important, the additional work (which is used in sending out sound 
waves) is all done before the boundary layer separates from the rear of the sphere. 
The work done on the fluid is increased twofold as a result of compressibility in the 
limiting case of infinite acceleration, but by a smaller factor for a finite acceleration. 
These results apply when the Mach number is low; there is some discussion of how 
they are affected by increase of Mach number. 


1. Introduction 

THE general effect of an incompressible, inviscid fluid on a sphere moving 
through it is to increase the apparent mass of the sphere, that is, in the 
dynamical equations of motion of the sphere a term of the form m-—+m’ 
occurs, where m is the actual mass of the sphere, and m’ is half the mass 
of the fluid it displaces (1). The quantity m+ m’ is called the virtual mass 
of the sphere. Further, because of the properties of an incompressible 
inviscid fluid, when once a sphere in such a fluid is given a velocity, it will 
continue to move with this velocity indefinitely. There will be no dissipa- 
tion of energy, this being caused, in general, by viscous forces which are here 
neglected. The motion of the fluid at any instant is independent of the 
previous motion of the sphere; because of this the sphere cannot transfer 
its energy to the fluid. The steady state of such a system is established 
immediately the sphere is set in motion, since disturbances in a perfect fluid 
are propagated instantaneously. 

For the motion of a sphere in a compressible, inviscid fluid different con- 
siderations apply. As viscosity is still supposed zero, there can still be no 
dissipation of energy by this means. However, because the fluid is now 
compressible, the problem becomes essentially one of acoustics. Distur- 
bances in this fluid will not be propagated instantaneously, but with the 
speed of sound, and hence whenever a sphere in such a fluid is set in motion 
in any manner, the steady state of the system will not arise immediately. 
but will take time to establish itself. Once the sphere has been set in motion, 


a continuous force must be applied to it to maintain its motion until the 


[Quart. Journ. Mech. and Applied Math., Vol .V, Pt. 1 (1952)] 
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~ 


steady state is established. If no such force is applied the effect of the fluid 








iE is gradually to reduce the sphere to rest. In this process the kinetic energy 
ff the sphere is transferred to the fluid, which acquires both kinetic and 
potential energy and radiates them to infinity. However, if the motion is 
maintained, the concept of virtual mass can be given a meaning as the total 
| work done to set the sphere in motion with uniform velocity, divided by 
half the square of the velocity of the sphere: it is this quantity which will 
he investigated here. 
Since through its motion the sphere causes acoustic energy to be radiated 
d te to infinity, in addition to building up the kinetic energy of the steady 
ti motion, the energy of the fluid, and hence the virtual mass of the sphere, 
ee may be expected to vary with the method used to impart to the sphere its 
ere, elocity. This is, in fact, found to be the case. Two methods of generating 
th the motion are considered. Firstly the sphere is given a certain velocity by 
“ae the application of a suitable impulse, and this velocity is then maintained 
onstant by the application of a suitable force. Secondly, we examine the 
se of the sphere accelerated gradually from rest to a certain velocity, and 
vhich is then maintained constant by applying a force. This method of 
ying generating the motion may be regarded as obtainable from the first by the 
the principle of superposition, since it consists of the limit of a series of small 
mi mpulses applied at short intervals of time. In both motions the pressure 
1ASs nd density are assumed to obey the adiabatic law, and in deriving 
1ass the equation satisfied by the velocity potential, the squares of any velocities 
‘ible have been neglected, since all velocities, and in particular the velocity U 
will f the sphere, are assumed to be considerably less than the velocity of 
ipa sound a. Hence M, the Mach number of the flow, is very much less than 
here nity. Points in space are referred to axes fixed at the centre of the sphere 
the ind moving with it. and velocities are measured relative to axes fixed in 
sfet space. Both the cases of motion of the sphere are discussed assuming that 
hed the velocity potential ¢ satisfies the acoustic equation. The energy of the 
luid fluid. and the work done on the sphere to maintain its motion, are calculated 
ndependently and (of course) found to be equal. The virtual mass is also 
on found, and in the motion generated impulsively the coefficient of m’ tends 
» no to twice its value for the case of a sphere moving in incompressible fluia 
now ith constant velocity. For a general acceleration law the expression for 
tur this coefficient is too complex to give any information, but the particular 
the 4 case of uniform acceleration is discussed. This shows that the ratio of the 
tio! oefficient of m’ in the compressible case to the coefficient of m’ in the 
ely ncompressible case, drops fairly rapidly from two to one, as the radius of 
ion the wavefront when the final velocity is reached increases relative to the 
the tadius of the sphere 
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A more accurate equation for ¢ is then considered and the weights of the 





terms in this equation are compared. This equation is essentially the | a 
acoustic equation with two extra terms added representing the convection oi 
of sound, one due to unsteadiness in the flow and the other due to the oe 
instantaneous flow pattern. If the variation in a, the velocity of sound in oa : 
the fluid at points where its density is p, is considered, then a can be — 
expressed in terms of ay, 0¢/ct, and qg, where a, is the velocity of sound adi 
through the fluid at rest (and of density py) and q is the particle velocity. 21 
This introduces two further terms into the equation, which now consists of ses 
the acoustic equation with speed of propagation a,, and four extra terms. . 
The flow is considered to be rapidly varying, for if this is not so the equation vans 
for ¢, correct to a first approximation, reduces to its form for an incom- aa 
pressible fluid. With this assumption regarding the flow, and remembering ng 
that all velocities are considerably less than ay, the weight of the term due ~ 
to unsteadiness in the flow compared with the term due to the instantaneous dom 
flow pattern is as 1: M. Also, of the two extra terms introduced by con- wie 
sidering the variation of a, the one involving é¢/ct is of the same weight as — 
the term due to unsteadiness of the flow, and the term involving ¢? is of = 
weight M compared with the same standard. Hence, to a first approxima- the 
tion the term involving q¢? and the one due to the instantaneous flow pattern _ 
may be neglected. The equation for ¢ now consists of the acoustic equation a 
with two terms added, and both of these are of weight M compared with pro] 
a weight of 1 for the terms of the acoustic equation. Using this 
modified equation, the form of 4 for both the impulsive motion of the sphere and 
and its gradual acceleration is calculated, but the complete solution is not fixe 
given as the analysis became too complex. However, the form of the energy T 
of the fluid is calculated and the extra terms in the expression for the energy 
due to the extra terms of order M in the equation for ¢ are found to be of 
order M?. The inclusion of the two further terms in the equation for 4, 
which are of order M?, would also give terms in the energy of order .M? (and whe 
higher orders). Hence the energy cannot be calculated to any greater 
accuracy by only including the terms of order M in the equation for 4, since The 
these give rise to no terms of order M in the expression for energy, but only con 
to terms of order M?. and more terms of this order occur due to the terms equ 
of order M? in the equation for ¢. 
In practice the flow due to the motion of the sphere would be irrotational 
only in its early stages, for boundary-layer separation will occur and then : 
eddies will form behind the sphere. For example, in the case of a circular whi 
cylinder which is uniformly accelerated (2), boundary-layer separation 
occurs when it has moved a distance of half its radius. (For the cylinder and 
subject to impulsively generated motion the distance moved is one-third of 
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the radius.) These limits indicate that for values of the acceleration for 
which the compressibility affects the virtual mass at low Mach number, the 
steady flow is reached before viscosity has any appreciable effect. Hence 
the expression found for the additional energy, which must be made avail- 
ible for radiation to infinity in cases of very rapid acceleration of the sphere, 


is always accurate whenever this energy is in fact significant. 


2. Impulsive generation of the motion of the sphere 

Consider a sphere of radius b, at rest in a compressible, inviscid fluid 
which is also at rest. The sphere is given an impulsive velocity U (assumed 
considerably less than the speed of sound a) in the direction Oz, and this 
velocity is then maintained constant by the application of a suitable force. 
The application of this force will result in a certain amount of work being 
done, and the effect of the force will be to maintain the energy of the sphere 
at its initial level after the impulse has been applied. The fluid will gain 
energy, both potential and kinetic, as a result of work being done on the 
system. Some of this energy is dissipated to infinity by the fluid. Let O be 
the centre of the sphere, and P a point fixed in space, with coordinates (r, 4) 
referred to O as origin. The motion is irrotational and therefore a velocity 
potential d exists which, to a first approximation, satisfies the equation of 


propagation of sound, 1 ed 
V* : irs (1) 
a* a 
ind which is such that the velocity of a fluid particle at P, referred to axes 
fixed in space, 18S q orad ¢. 
The boundary condition on the sphere is 


/ 


CO 7 
U cos @, (2) 
a) 
when ? b. for all values of 6. The initial conditions of the motion are 
hd Cd ot 0 when t 0. 


The motion is symmetrical about Oz and so (guided by the boundary 
condition (2)) we assume a solution d Reos@. where R R(r,t). Then 
equation (1) becomes 

eR  2e0R 2R 1 eR 


) 2 9 9° 
or yr co ' a‘ at* 


(3) 


) 


while the boundary condition (2) is now 
CR/c) U for r b, 
und the initial conditions on R are 


R OR /ct QO at tft 0. 
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Hence, writing (3) in the notation of Heaviside’s operational calculus, it 


becomes @R 2dR 2R- p? 
— ine eee - R. 
dr? ° r dr i a 
Substituting R = Yr-}, this latter equation reduces to 
ay 1dY AT a y—0. 
dé rdr \4e ( @ 


Bessel’s equation for functions of imaginary argument. The two indepen- 
dent solutions of this equation are 

Y = K,j(pr/a) and Y = J,(pr/a). 
The appropriate solution for an outgoing wave is K;(pr/a), since in the 
asymptotic expansion of this function there occurs a term e~?"@ indicating 
an outgoing wave. Hence, for some A independent of r, using the known 


expression for A;(pr,a) as an elementary function, 


ova y 
R at) af =a) pria (4) 
» , 2 
<p ? | ia 


rom the boundary condition on the sphere the constant A is evaluated 
and the result is 


Uab? ‘a 
R ab? p ( a } pr—D)la_ (5) 


(p*b?+-2apb+2a*)\r — pr? 


Interpreting this operational form, we find 


Ub? cos 6 


dy2 


ar" 


a 
| l e—alt-(r—b)/a "C08 5 lt (r b) a|-- 
) 


l . @ 3 
(2r—b)e—all--Hallbsin —[t—(r—b)/a]|, (6) 
b b : 
for ft (r—b)/a. that is for r < b+at. d 0, for t < (r—b)/a, that is for 
) b--at. It follows that 

Ch Uh? cos @ 


Cc? ry? ) 





l = ; 
e—"” eos uw -t 5 (or bye ’ sin w 4 


Ubeos@ ? 


Ss le ’ eos w—e- sin w|, 


m 
where w = alt—(r—b)/a}/b. 


It is easily verified that the above value of ¢ satisfies the boundary con- 


dition on the sphere. Further, when r = b-+-at then w = 0, and hence 
od Ub cos 6 
or =—ssb++-att 


This is the value of the normal component of velocity at the wavefront, 
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calculated from the value of ¢ for the motion inside the wavefront. Outside 
the wavefront the fluid is at rest and hence ¢¢/er = 0. Thus there is 
, discontinuity of the above amount in the normal velocity at the wavefront 
that is essentially a shock. This discontinuity decreases as ¢ increases, 
vecording to the ordinary laws of spherical attenuation. However, this 
effect will not be exactly correct for any length of time as the neglected 
squares of velocities will begin to have an accumulative effect. The tangen- 
tial velocity at the wavefront is continuous, being zero at both sides of it. 
[t is also easily verified that the above value of ¢ satisfies the boundary 
condition at the wavefront namely, a é¢/ér od / ct. 

The disturbance terms in the expression for ¢, that is the departure of ¢ 
from the final steady value — U’b* cos 6/2r?, are rapidly damped for large t. 


\t any fixed point (7, @) they fall exponentially with f. 


3. Energy of the fluid, calculated both from the work done by the 
sphere and from the potential field 


From Bernoulli's equation, neglecting squares of velocities, 
P Po Po Ch ct. 


where p is the pressure at a point P in the fluid, and p, is the pressure at 
points where the fluid is at rest 

Hence the thrust, X(¢), on the sphere in the negative Oz direction is given 
by - 


X(t) | | p|,-, Cos 6 27h? sin 6 dé 


27 py 6° | cos 6 sin 6 dé 
< r=b 


ot 
0 ; 
Aor, ab? Le a/b eos ath. (7) 
The initial drag coefficient is 
X(0) 8 
Lp, U2? 3M’ 


Po 


where M is the Mach number. The drag is periodic, because of the presence 
of the term cos at/b, and is also damped exponentially with time. Thus it 
dies away and the steady state of the motion is then approached, when the 
sphere continues to move with constant velocity and is subject to no drag. 

To maintain the velocity of the sphere constant, a force must be applied 
to overcome the drag. If m 3771p, b® (= half the mass of the fluid 
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displaced by the sphere) the magnitude of the force which must be applied 
is given by F(t), where 


2m'a 


F(t) 


Ue-“! cos at/b. 
b ; 


and the force acts in the positive Oz direction. 
The work done by this force in time ¢ is given by Q, where 


t t 
Q(t) | F (7) dz | F(z7)U dr 

T=0 0 

t 
uit h 
“MA +75 e-27lb eosar/b dr 
“el 
0 
m'U?| 1 +-e-@ sin at/b—e-" cos at/b|. (8) 


“quation (8) will now be obtained by an alternative method. 
The kinetic energy 7 of the fluid is given by 


T , | py dv, 
Vv 
where J’ is the volume occupied by the fluid and dv is an element of this 
volume. If p = p,(1+-s), where s, the condensation of the fluid, is small, 


then . 
T = 4py { gidv. 
Vv 


At time ¢ the volume of fluid which possesses any energy lies between the 


sphere and the wavefront, that is, between r band r bat. 
Therefore : 
[= Po | q* dv 
V 
b+at am 
spy | | q?2rr? sin 6 dédr. 
b 0 


The value of ¢ is of the form R(r.t)cos 8. where Ris known. Hence 
y* (eo R/er)? cos?6-+- R? sin26/r? 


b+at 


. 9 2 2 
sy 2a? () | - | dr. (9) 
3\ Or 3 r 


b 


and therefore aT 


after integrating with respect to 6. 
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The boundary of the volume V enclosing the fluid which has a motion, 
and hence kinetic energy, varies with time, and thus in calculating d7'/dt 


there will be a term due to this effect, 


dT df 
Po (orad 4)? dv 
dt "ae 
. 
op 4R2 22 /@R\? 
Mo grad d. grad — dv+-mp,ya|——-+ 
. ot 3 3 \ er colieas 
J 
: Ch = a 
™ | grad }. grad — dv+ §apyab?U®. (10) 
ot 3 


« 


} 
The term in d7'/dt due to the boundary enclosing the moving fluid changing 
with time is easily calculated from the form of 7' given in equation (9). 
The internal, or potential, energy W of the fluid, which is due to the 
compressions or rarefactions occurring, is given by 
, ee , * 1 (ed\? 
i 2 Po a*s* dv > Po | re ape dv. 
: J a*\at 
V 
provided that the initial and final volumes of the disturbed mass of gas are 
the same. Hence 


1 (ed\? 
YW > Po (*) 2rr? sin 6 dédr. 


a*\ ct 


Sut cd ct is of the form F(r,t)cos 6, and hence 


ft 


. 


Vi =Pe | 12? cos*# sin @ dédr 
a= 


<7 Po F2r2 dr. (11) 
It then follows that 


( » ad ff [ed\? 
Ai Po | ) dv 
dt 2a? dt d ot 


moL 294 5) 

) CO C-@D «TT —. 19 6 
Po [ob 26 ay 4 2mPOr papal 
a“ J ot at 3a 

J 
) y Ch ( “4 
Po - — dv-+ 3p, ab?U?, (12) 
a" ct ct? sa 

} 


the expression for W from equation (11) being used to calculate the term 
indW dt due to the boundary of V changing with time. 





=o 
la 


A. L. LONGHORN 


Adding equations (10) and (12) then 
d 
dt 


a ot ct? 


ay ' Cd 1 ede aes 
(T+W) = py | (arad g.grad :) dv+imp,ab?l™. 
C 
I 
By Gauss’s divergence theorem. 


. 


f . d we 
grad d. grad “ dv Op Op ds | od V2dh dv. 

: ot J ot on J a 

V S V 


where S is the surface bounding the volume V, n is the unit inward normal 
to this surface, and ¢/én denotes differentiation along this normal. 


Applying this result, and using also the equation for ¢, we then have 


d __ > 4 Ch Cd - erro 
(7'3- W) Po dS Sarp, abl 2. 
dt J ot on a 
S 
The surface S consists of the spheres r bh and r b+at. Onr=b, 
Od /en Cd /er Ucos@, éd/et Uacos 0e-“ cosat/b. On r hat, 
Cd /en 0d /er Ubecos@/r. ed/et acd /er Uab cos @'r. 
Hence d 


(7+W) Sap, ab?Ue—-“ cos at/b. 
dt 
and therefore the total energy of the fluid at time ¢ is 77+ W. where 
t 
T+W tarp, ab? [ a/b cosat/b dr 
0 
atlb 


arp, ab?U>| 1+ sin at/b—e-“” cos at/b| 


m' LU | -+-e-“ gin at/b—e-@! cos at b|. (13) 


From equations (8) and (13) the values of the energy of the fluid. calculated 
by the two methods, are seen to agree. 

Initially, when the fluid is at rest and the velocity U’ has been imparted 
to the sphere. the total energy of the system is }mU?. In the subsequent 
motion the force applied to the sphere maintains its velocity, and hence its 
kinetic energy. constant, but because of the work done by this force the 
energy of the system increases, that is the fluid acquires energy. By the 
calculations from the equations (8) and (13), the work done by the force is 
found to be the same as the energy acquired by the fluid. In an incom- 
pressible inviscid fluid, when the motion is set up by the application of an 
impulse to the sphere, the liquid acquires a motion immediately and a 
steady state exists. The motion of the fluid is dependent on the motion of 
the sphere and ceases with it, and further, the fluid offers no resistance to 


the motion of the sphere. Thus the initial impulse applied will give the 


sphere its kinetic energy 4mU? and also give the fluid its kinetic energy 
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sphere, and since there is no dissipation of energy, the energy of the system 


No further force will be required to maintain the motion of the 


will remain constant. Since the fluid is incompressible it can possess no 
potential energy as compressions cannot occur. On the other hand, equa- 
tion (8) or (13) shows that, when compressibility is taken into account the 
total energy ultimately imparted to the fluid is m’U?, namely twice as much 


is when the compressibility is ignored. 


4. The sphere subject to a finite acceleration acting for a finite 

time interval 

\ sphere of radius b moves with velocity U(t) through a compressible, 
nviscid fluid. For all values of the time t, U(t) is considerably less than the 
speed of sound a, and when ¢ is zero the sphere and the fluid are at rest. 
Let O be the centre of the sphere and P a point fixed in space, with co- 
dinates (7, @) referred to O as origin. The sphere moves in the direction Oz. 
‘he motion is irrotational and therefore a velocity potential ¢ exists, and if 


the velocity of a fluid particle relative to axes fixed in space be defined as 


q = grad 4, then to a first approximation ¢ satisfies the equation 
; 1 od 
V2h ee te (14) 
a* at 
[he boundary condition on the sphere is 
U cos (15) 


when ? b. for all values of 6. 


Cd ct 0 when ft 0. 


Because of the form of the boundary condition (15), we assume a solution 


The initial conditions of the motion are ¢ 


of the form db R(r,t)cos @. 


Then, since the motion is symmetrical about Oz, equation (14) becomes, on 
substituting for ¢ and in the notation of Heaviside’s operational calculus, 
dR 2dkR 2R ? 
: r= ak. 
dr? } d) = a- 
Solving this equation and choosing the appropriate solution for outgoing 


waves from the sphere we find 


7a\4/1 a 
2p , 


where B is independent ot r. 
Using the boundary condition ¢R/ér U when r b, then B is deter- 


mined, and Vab3 ra 
ao Pp T 


R 


p(r—b)ja 


p*b?+-2apb+ 2 
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Interpreting this by using the Faltung integral result, we have 


t 
l cree ; 
d ; ad | | (t —r)e alr—(r—b)/a}/b cos al - (r b) a| h dr | 

r 

0 
t 
L fas ce 
4+ (r b) cos 6 | U (t—r)e—ar-" al sin al r—(r—b)/a\/b dr, (16) 

r2 


0 
if (r—b)/a < t, that is if r << b+at, and d 0 if (r—b)/a > t, that is if 
r > b-+at. Changing the integration variable, we have the alternative 
form of the solution 


t—(r—Db)/a 
» 


ab | 


db a= OGRE U (r)e—all-7r-Yialb cos alt— 7—(r—b)/a]/b dr 


. 
0 
t—(r—b)/a 
. 


l : ; : 
(r b) cos 0 U(r)e—alt-7-¢ Dallb sin alt—r (r—b)/a|/b dr (17) 


ifr < b+-at and d = Oifr > b+at. 

It is easily verified that this value of ¢ satisfies the boundary condition 
on the sphere, and also the boundary condition at the wavefront, namely 
a od/er é¢/ct. Further, at the wavefront r = b--at, éd/er = 0, and 
hence there is no discontinuity in the normal velocity there. 

The solution for this case of finite acceleration of the sphere could have 
been obtained from the solution for the impulsive generation of the motion 
by a principle of superposition. The gradual acceleration for a finite time 
interval can be considered as the limit of a series of small impulses, applied 
at successive small intervals of time. This gives the result (17), and further 
it is easily verified that if U(t) is a constant. then this solution reduces to 
the solution for the impulsive motion. 


5. Calculation of the energy of the fluid, both from the work done 
by the sphere and from the potential field 
From Bernoulli’s equation, neglecting squares of velocities, the pressure p 


at a point P is given by DP = Do— Py od /ét. 


Therefore X(t), the thrust on the sphere at time ¢, is given by 


( [ed a 
X@ = — | of | cos 6 2mb2 sin 6 dé 
« ot r=b 
0 
t 
“aU 
tarp, ab* | on email 7 cos a(t—r)/b dr. (18) 
dr 
0 


and it acts in the direction of negative Oz. 
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To maintain the velocity of the sphere a force must be applied which will 
neutralize the drag. This will have the same magnitude as the drag but act 
in the opposite, that is positive Oz, direction. The work done by this force 


in time ¢ will be given by Q(t), where 


t t 
Q(t) ( X(t,) dz ( X(t,)U(t,) dt, 
t;=0 0 
: ‘ 
br p, ab? } U(t,) iP “h—7) cos a(t,—7)/b drdt, 
i 0 
 » Dan 
= 2 | ol “h—™eos fa(t,—7)/b} U(t,) drdt,. (19) 
0 6 


Equation (19) will now be obtained by an alternative method. The kinetic 
energy, 7’, of the fluid is given by 
Zz $ pq” dv = Po q* dv po | (grad d)? dv. 
i Vv 
where V is the volume occupied by the fluid. 
But R(r,t)cos@ and therefore gq? = (@R/ér)?cos*@+ R? sin?6/r?. 





Thus , >\ 2 2 
a Ca 4. 
1 S09 | cos? -_ — sin?6| dv. 
is 2 Cr | ia 
J 
\t time ¢, the volume of the fluid which is in motion lies between the sphere 
b and the wavefront r = b+-at. Thus 
b+at am 
” , OF CR ws 2c 
| 1 09 | ( cos?@ + — sin?6 | 27r? sin 6 ddr 
cr ie 
b 0 ‘ 
b+at 
f [2r2/eR\?  4R? 
T Do ) -| dr. (20) 
j oa \o 3 
b 


The volume occupied by the fluid in motion varies with time, and thus 
in calculating d7'/dt this variation must be considered. The form of 7’ from 
equation (20) will give this result conveniently. For 

iT , a 


15 (orad 4)? dv 
dt or ot | 


"a 


F C db ‘2r2 oR > E 4R? 
Do grad ¢. grad — dv+-mp,a — mi ‘ 
t r=b+alt 


3 \ér 3 


« 


J 
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But ¢d/er 0 when r = b-+-at and hence @R/ér = 0 there. Also d = 0. On? 
and hence R 0. for r b+-at. The second term above is therefore zero 
and hence 
dT a od 
- Po | grad d. grad — dv. (21) 
dt H ot 
"g Hence 
The potential energy of the fluid due to compression and rarefaction is W, 
where * Log\2 
y C@m\~ 
i Po dv. ai 
2a? ot [here 
. 
But ¢d/ét (eR/ct)cos 8, and hence 
b+at 
: o, { (eR\? ) f OR\? 6 
Vi fe. cos?6 dv ~ | ( cos?6 2rr? sin 6 déd) 
2a* J \ at 2a* . ot 
V b 0 
b+at 
‘ : ‘ - 
27 ,(cR\2 Thus 
: Po isi dr. (22) 
0" ot in (1¢ 
main 
Using (22) to calculate the term in dW’/dt due to the fact that V changes | just 
with time, we find (‘5 
1W ] > d 2 . A 2d, 9 R 2 at til 
¢ Po § (: de a | OP OP gy 1 =7Po | ya C# 
dt  2a®dt J) \ét a2} et ét 3a | \ at) |ecpeat er 
Vr V wher 
But 0¢/et = 0, and hence @R/ct = 0, when r = b+at. Thus the second 
term above is zero, and 
dw Po } od ¢ 2h / (23) [his 
: - _ av. . Zo ° 
dt a2 J] at é {2 inte: 
V thar 
Adding equations (21) and (23), 
pen ’ od 1 he 
(7+W) = py grad ¢d. grad | dv 
dt 4 ae ada oF oe 
V This 
ig od od : inte 
ei oe, dU 
ot on 
Ss bou 
where S is the surface bounding the volume V, n the unit inward normal stat 
to this surface, and ¢/én denotes differentiation along this normal. The I 
above follows on applying Gauss’s divergence theorem and using the 
equation satisfied by ¢. In this case the surface S consists of the spheres Q( 
r= bandr = bat. 
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On) bat. é6d/ét 0. On? b. ed/en od/cr U(t)cos 6 and 
Cb F dl 
acos@ e—al 7b cos alt T) bdr. 
ct ar 
Hence 
f 
i tor “dU 
(7 W) Po abl (t) | oe enaut 7b eos a(t—r)/b dr. 
dt 3 ; dr 
0 


Therefore the total energy of the fluid at time ¢ is 7+ W, where 


we . . ° dU 
fl HI $77 py ab? ( (0; = i cos a(t; T) b drdt, 
0 . 
2m'a [ f dl _ 
r | U(t, = ah—7 eos a(t, —7)/b drdt,. (24) 


[hus this value of the total energy of the fluid agrees with the value found 
n (19) for the work done on the sphere to overcome the drag and hence to 
naintain its velocity. To generate the motion of the sphere more work 
nust be done, of amount ml?/2, where m is the mass of the sphere. 
Consider a sphere which is accelerated in an arbitrary manner from rest 
it time ¢ = 0 to a velocity LU (t) l), at the time ¢t = f,, after which this 
elocity is maintained constant. The drag on the sphere at time ¢ is D(t), 


VW here 


Dit) —_ e—at—r)1b eos a(t—7)/b dr. (25) 


Chis is true for all t, but if ¢ is greater than t,, then the upper limit of the 


integral is t, since dU’ /dr, and hence the integrand, vanishes for 7 greater 


than ¢ 


0° 


Hence for all ¢ greater than fp, 


2m'a f dl 
Dit) — e-al—7)l5 eos a(t—r)/b dr. (26) 


[his expression contains the factor e-“”, and if this is removed from the 
integrand then the latter is bounded in the range of integration provided 
IU dz is bounded. Hence. because of the factor e-“”, assuming dU /dr is 
bounded, D clearly tends to zero as t increases. When this occurs the steady 
state of the motion is approached. 

The work done in time t is given by Q(t), where 


2m a 


b 
ny 
0 


dU 
€ 


Q(t) D(t,)U (t,) dt, ; 
at 
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0 


aly “7 cos a(t, —T) b drdt,. 
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jf ow 
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This expression for general values of dU /dr is too arbitrary to yield much 
information, and hence a special case will be considered, namely a linear 
acceleration of the sphere. 

Suppose the sphere is accelerated according to the law u = At. where dis 
a constant, from ¢t = 0 to t= ¢, and then u = uw, = At, is maintained 
constant. Substituting these values in the general expression for Q(t), the 
work done in time ¢, then 


a & 
2m'a f 
(V(t) j | | r*t, € at 7 cos a(t, t)/b drdt, : 
} 

0 0 

tt 

Smal, [ £ 

; = | | Ae “4-7 Cos a(t, —7)/b drdt,. 
iy 0 


Performing the integrations we find 


m OE b2 


V(t ) ( pe —alolb 


cos at, /b—e-“e” sin at,/b) + 
5) 


z att 


2h 


! (e-4¢—0) sin a(t—t,)/b—e-@ sin at/b)|. (27) 
ato 


When ¢ is large compared with ¢,, then 
m'U2 


») 


Q(t) 


b2 
1 + — (1—e-“ cos at,/b—e-“o! sin at,/b) (28) 
at) 

nearly. 

Hence the virtual mass of the sphere when ¢ is large compared with f, is 

m+-m'| 1 +-(b?/a?#2)(1—e-“!” cos aty/b—e-“o” sin aty/b) |. 

In the case of a sphere moving in an incompressible fluid the virtual mass 
is m+m’. Hence if the ratio of the coefficient of m’ in the expression for 
the virtual mass for uniform acceleration in compressible flow to the coeffi- 
cient in the expression for virtual mass in incompressible flow be y, and 
at,/b be x, then 


y 1 4 — (1—e-* cos x—e-“sin 2). (29) 


Considering the graph of this function, we find that as x > 00, y — 1, this 
being the case of a very gradual acceleration over a long interval of time. 
As x > 0, y > 2, and this represents the case of a steep acceleration for a 
short interval of time. In the limit this gives the impulsive case. In all these 
cases u, is regarded as constant and the acceleration A and the time-interval 
t, are varied to give quicker or slower accelerations. 


In practice, boundary-layer separation will occur and then eddies will 


begin to form behind the sphere. Some calculated results on the distance 
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moved by a circular cylinder of radius b before separation occurs are avail- 
‘ble (Goldstein and Rosenhead (2)). When the motion of the cylinder is 
venerated impulsively separation of the boundary layer begins at its rear 
stagnation point after it has moved a distance of 0-32b. For uniformly 
.ecelerated motion of the cylinder separation again begins at the rear 
stagnation point after the cylinder has moved a distance of 0-526. Similar 


results may be expected for the sphere. Therefore the parameter x, in 





Fic. 1. 


terms of which the results are expressed in equation (29) and Fig. 1, can 
reach a value about }./~-', which is fairly large, before the accuracy of the 


solution is at all affected by the action of viscosity. 


6. More accurate equation for the velocity potential 
When squares of velocities are not neglected the equation satisfied by ¢ is 


avd — 41°! 1 orad dd. grad 49. (30) 
ot ot - 
Further, the velocity of sound a varies with the density of the fluid, and if 


the latter behaves according to the adiabatic law then 


’ Ch 
a= a- ) 1) Hy 1)q?. (31) 


ct 


vhere a, is the velocity of sound through fluid at rest and of density pp. 
In previous work, the acoustic equation has been used as an approxima- 
tion, the last two terms of (30) and the variation in a being neglected. These 
extra terms will now be included and their relative weights considered for 
the case of certain flows. Consider a sphere of radius b, moving with velocity 


considerably less than a). through a compressible, inviscid fluid. The 
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velocity of a fluid particle relative to axes fixed in space is q, where q is | 


grad ¢, and q is the magnitude of q; ¢ is also considerably less than ao, that 
is, the Mach number M (= q/a,) of the flow is conside rably less than one, 
From (30) and (31) the equation for 4 which is to be considered is 
9 9 ‘ 
‘i on Oh . oq” . 0b —. sain 
ap V°dh —— - } orad d. grad 4g? (y Ns V2¢6-4 iy 1)qg? V7d. 
C 


/ 
at C 


Now o/éx = 1/b, éd/éx = q and hence ¢ = bq. If also G/ct = n then, 
az V*d = az q/b, od ct? = n2qb, 
0q7/et = nq", grad 4. grad(q?/2) = q3/b, 
(y—1)(6d/et) V*d = ng? and (y—1)(q?/2)V2d = q3/b. 
Hence the relative weights of these six terms are 
1, n*b?/as. nbg/az, q?/az, nbq/a2, and q?/a?. 
or 1, n*b?/a%, Mnba,, M*, Mnb/a,, and M2. 
If nb/a, = 1. then the relative weights of the terms are 
1, 1, M, M?, M, and M?, 
Since M is considerably less than one the terms of weight M?* are considerably 
less than those of weight WM. In the acoustic equation terms of order M and 
M* are neglected compared with those of order one. However, to obtain a 
more accurate equation for ¢, the terms of order M are included. 
Hence the equation for 4, to an approximation which neglects V2 but 
not M, becomes 
l od 1 cq? , (y—1) ed 


az ct? ~~ az et az at 


V4 V4. (33) 


~~ 


[It has been assumed above that nb/a, = 1, that is that nb = a,. Now 
o/et = n and hence 
(wo/ox)/(0/et) & q nb = JA M l. 

3ut (we/cx)/(e/et) is a measure of the rate of variation of the flow. and for 
rapidly varying flow this ratio is considerably less than one. This is just 
the result that follows from the assumption nb = ay, and hence the weights 
of the terms as here calculated hold for the case of rapidly varying flow. 

For slowly varying flow g/nb = 1 or q/nb l. Ifthen q < ap, it follows 
that a)/nb > 1, that is nb/ay <1. With this result the equation for ¢, to 
an approximation neglecting second-order terms, becomes V2é = 0, the 
equation for incompressible flow. 

An attempt was made to solve the equation (33) in the case when the 
sphere was given an impulsive velocity U, which was then maintained 
constant. The value of ¢ obtained for the same problem when ¢ was 








assum 


oq” ct 


where 
functi 


Eq 


where 
of arg 
when 
the fl 
Lege! 
to inc 
Th 
from 
will b 
equat 
in th 
bette 
equat 
term 
terms 
while 
sectic 
Th 


encol 


509 







































THE UNSTEADY, SUBSONIC MOTION OF A SPHERE 8) 








418 | . sumed to satisfy the acoustic equation was used to calculate the terms 
that aq? /et and (é¢/et) V*¢d in (33). The resulting equation is then of the form 
me. 
1 od ’ 7 
V*¢ =e (7, t).Py(cos 6) + F,(r, t)P,(cos 6), (34) 
0 
. where F,(r, t), F,(r, t) are known functions of r and t, and P,, P, are Legendre 
(32 functions. 
Equation (34) has a solution of the form 
@= MA, P,+A4,Ah4+MA,P,, (35) 
where Ay, A,, A, are functions ofr and t, and P,, P,, P, are Legendre functions 
} ofargument cos @. The solutiond = A, P,(cos @) is that which was obtained 
when it was assumed that ¢ satisfied the acoustic equation. The energy of 
the fluid can be calculated and because of the orthogonal properties of the 
Legendre functions the extra terms in the expression for the energy, due 
to including the extra terms in the equation for ¢, are of order M?. 

The same form of solution will hold for the sphere gradually accelerated 
from rest, and hence again the additional terms in the energy of the fluid 
will be of order M?. However, if the further two terms are included in the 

ably equation for ¢, then still more terms of order M?, and higher orders, occur 

and } in the expression for the energy. Thus from energy considerations, if a 

in a better approximation than the acoustic equation is to be used, then 
equation (32) is required, as the terms in it of order M? contribute a 

but term to the energy of the same order as that contributed by the two final 
terms in equation (33). The labour required to do this will hardly be worth 

i while: it is at least encouraging to note that the virtual masses obtained in 

i sections 3 and 5 are correct if terms of order M? are neglected. 

Vow The author wishes to thank Professor M. J. Lighthill for his help and 
encouragement during the investigation of this problem. 
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THE DIRECT USE OF GREEN’S METHOD FOR 
SUPERSONIC POTENTIALS 
By J. OKEEFFE (Department of Mathematics, King’s College, London) 


[Received 13 February 1951] 


SUMMARY 


This paper gives an investigation of Green’s method for linearized steady super- 


sonic motion past a thin wing without making any appeal to Hadamard’s theory of 


the finite part of an improper integral. 


1. Introduction 

In the classical hydrodynamics of an incompressible fluid the velocity 
potential can be expressed by Green’s formula as equivalent to that due 
to a suitable surface distribution of sources and sinks. 

The corresponding problem for a compressible fluid admits a very similar 
solution if we restrict ourselves to small steady disturbances superimposed 
on a steady uniform flow, provided that this basic flow is subsonic. 

If, however, the basic flow becomes supersonic, the character of the 
problem is scnighaely changed; the partial differential equation for the 
disturbance potential ¢ becomes hyperbolic instead of elliptic; and, as a 
consequence, various special and artificial devices are employed to con- 
struct the analogue of Green’s formula. Of these, the most famous are the 
methods devised by Volterra and Hadamard and described by the latter 
in (1). 

The purpose of this paper is to show that a simple and direct solution 
of the problem of steady supersonic motion over a thin wing at zero 
incidence can be constructed by methods strictly analogous to the classical 
methods of Green. The advantage of such a solution is that it avoids any 
appeal to Hadamard’s general theory of the ‘finite part’ of an improper 
integral by showing explicitly that the apparent infinities exactly balance 
one another. 


2. The analogue of Green’s method 


Green’s method for expressing classical potentials by means of surface 
integrals consists in applying the formula 


IS bi)-S) mle st)-23) ea) -F3)] 


en 


2 | | | {sa(- ait fi (#)} dgana 


(Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 1 (1952)] 
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where the volume integral is taken over any closed regular region D 
within which the integrand is regular, the surface integral is taken over 
the set of surfaces S which bound that region; 1, m, are the direction 
cosines of the outward-drawn normal to the surface element d=, 
r? = (E—2)-+(q—y)*+(E—2)%, 
2 g g 

and A + +>: 

o€ a 
to the region D, bounded externally by the sphere r? = R? (R > 1), and 
internally by the sphere r? = e? (e 1) and a given surface on which ¢ and 
its first partial derivatives are known. If ¢ is harmonic, the volume integral 





an. ? 
on” 


vanishes, the integrand vanishing everywhere in D,; the part of the surface 
integral due to the sphere of radius R is O(R-!) as R tends to infinity; 
and the part due to the sphere of radius « tends to 47 ¢(x, y,z) as « tends 
to zero; so that proceeding to these limits we obtain a formula for ¢ at 
ny point (x,y,z) outside the given surface in terms of the known values 
of ¢ and its first partial derivatives on that surface. 


The analogous formula for supersonic flow is 


[meta ef)-22)-aoat)-2a 


ener 


[jf e+()—H0) ae 


where the integrals are taken over any region D’ and its bounding surfaces 
S’ where the conditions of regularity hold, 


a = B-4(E—2)*—(q—y)*—(C—2), 


9 a2 a9 
ia Oo” - 


and L B? 


a yt a 
We now apply this formula to the problem of a thin supersonic wing at 
zero incidence in the plane € = 0, with the line 
¢ = @, € = X(A), 7 = Y(A) 
is leading edge, considering the two closed regions bounded by the hyper- 
> X, and the 


boloid s? = e?, the part of the plane € = 0 for which é > 
control surface’, i.e. the envelope of the cones 
B-*(é—X)?—(n—Y P-L? = 0. 

In the figure the closed part cut off on the upper part ({ > 0) of the 
ontrol surface by the hyperboloid s? = e? and the plane { = 0 is denoted 
by S}; the closed part cut off on the lower part (¢ < 0) of the control 
surface by the same hyperboloid and plane by Sj; the part of the 


sheet of the hyperboloid s? = «2, £—x < 0 cut off by the upper part of 
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the control surface and the plane ¢ = 0 by Sj; the part of the same sheet 
of the hyperboloid cut off by the lower part of the control surface and the 
plane € = 0 by Sz; and the closed parts of the upper and lower surfaces 
of the plane € = 0 bounded by the leading edge and the curve of inter- 
section of the sheet of the hyperboloid s? = e? (é—a < 0) and the plane 
¢ = 0 by Sy and Sg respectively. The leading edge in the figure is a 
straight line. 


la 
e \ 
4 
Fi / 
/ / ++ 
C # S ‘ 
/ ' ; 
\ ye 
Mg Pall 
c : ee a 4 
+ -* <i 
, / 
ale 4 ‘, 
/ - PA 
/ - ae, 
ae nose: ae a 
’ p <a TT L/S 2 
a hing Ss Ps 
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Since ¢ is a supersonic potential, L(¢) = 0, the volume integral vanishes, 
the integrand vanishing at every point in the two regions under considera- 
tion. Hence 


— re l l ed) { C (| l od) Fa: ] ] od) 
B l e< —_—-—- - t ») — 7 . - oe 
| | | | a.) 8 0&] : \¢ On ) s én} “? as) s elf 


the integral being taken over S; +S +Sj, or over S 





(1) 


1+S8,+Sy;. The 
supersonic data are the values of 4 and of its first partial derivatives on 
the wing, and the fact that d 0 on the control surface. 


3. Integrals over S; and S, 


Since ¢ 0 over Sf and S;, 


| | {Bre (")—m ohn 2 (4\\ as = 0 (2) 
Fier! \ c&\s On\s el\s){ 
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lerivatives of ¢ in directions parallel 


to the tangent planes to these surfaces must vanish. But at the point (€, 7, C) 


of the control surface 


und 


so that 


(l, m, n), oc (B-7(E 


the direc 


‘tion (6. m, n 


that point. Consequently 


l CO 
9 + 

| | - Vit 
Ss oc 


nN 
rs 
— 
~ 
= 
| 
~ 
us 


) is parallel to the tangent plane at 


C7) CG 


r 
op 4 ad = 0. (3) 


the integral being taken over the surface Sy; or over the surface Sj. 


4, Integrals over S; and S, 


Let 


i 


Bs cosh 6. ” Yy § sink 


1\@siny, and (—z = ssinh@cos y; 


then the equation of the hyperboloid in these coordinates is s e. The 


direction cosines of the outward normal to the hyperboloid at the point 


) 


/ 


4 


W here 


are 


Mt, Vl 


BoA ecoshé, A ‘sinh @sin x, A ' sinh @ cos x, 


A (B-* cosh?@-+-sinh?6)!, 


and the element of area on the hyperboloid is 


d> o*e* A sinh 6 dédy ; 


so that the surface integrals are transformed into the double integrals 


7 


ve 


B cosh 6 


| Cd d cosh 6 A 
)-+sin 


€ cé Be? 


sinh 6 cos 


taken between suitable limits Now 


ind 


so that 


f 
cc 
Be cosh 0 


Be sinh 6 
oo 
¢ . 
ecOsno sin y 
oo X 
ara 
€ cosh § GOS vy 
rales X 


e sinh @sin 4 


cc c 


db . : Oc 
-+-esinh@ cos y — = 


ih sin x 


€ ©7) ie 


l od ésinh @ cos ) 


) 


e- 


lcd dsinh@sin y ‘ 
x( cp @ in u 2%)! Bet sinh @ dOiy 


uw 


€ ¢€ 


Be(cosh 6—e-*), 
e(sinh @+-e~®)sin y, 


«(sinh 6-+-e-*)cos x: 


1] 
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where F=B = sin x na + cos x we 
a€ en ol 


thus our double integral can be written 


ep 








(3 +p+e-% F) Bsinh 6 dédy, 
B 


taken between suitable limits. 
Denote the values of 


| | e-% F Bsinh 6 dédy 


taken over Sj, Sy, by J, Jz respectively. 


Ti > a 
_ Iz |, Ig | < «Bmax F | | e~* sinh 6 dédy 
4 
27eB max F | e~* sinh 6 dé (y = cosh-(x/e)). 
0 
Thus I< 1s O(c log(1/e)). (4) 


Before considering the remaining parts of the integrals over Sj, Sy 
we first note that over any area of these surfaces 
i B ry 
| | (++ <5) sinh deat = | |\dsinh a — | ge? dd} dy. 
x 64 
where f, «, n(x), 94(x), are the limits for 6 and y corresponding to that 
area. 
We now break up the range of the variables 6 and y in the double integral 
over Sz into three parts as follows: 
(i) 6 = 0 to 0 = 6,, and xy = 0 to x = 27, where 


6, = cosh-!(k/e!) 
and k is any finite constant ; 

(ii) 6 = 6, to @ = 6, and xy = ato x = b, where a and b are limits for y 
depending on z and « and @,(y) is the value of 6 on the curve where the 
hyperboloid s? = e? cuts the upper part of the control surface; and 

(ili) 6 = 0, to @ = @, and x = b to xy = a, where 0,(x) is the value of 4 
on the curve where the hyperboloid s? = ¢? cuts the plane ¢ = 0, ice. 

6, = sinh-!(—zsec y/e). 
On Sz the corresponding range of integration is 

(iv) 6 = 6; to @ = 0, and xy = b to x = a, where 8,(x) is the value of 6 
on the curve where the hyperboloid cuts the lower part of the control 
surface. 


The range of integration (i) corresponds to a small cap near the vertex 
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of the hyperboloid S, when « is small and —z, n—y, ¢—z are all of order 
ein this range, so that 
4(0)—4(0) Oleh) (0O< 06 <4#,), 


where $(0) means ¢(a— Becosh6, y+esinh@sin xy, z+esinh@cos x) re- 
garded as a function of 6. Thus 

0, 0; 

| de-9 do = | {4(0)—¢(0)}e-6 d0-+-4(0)|1—e-9|" — Ket +-4(0), 


0 0 


so that the contribution to the integral over Sf from 
ei Op ° 
| | (<4 +¢)Bsinh 0 dédy. 
ol 


taken over this range, is 
27 Bd (x, y, z) + (8, )sinh 6, dx | Kel, (5) 
0 
In the range of integration (1i) 
(8) O(1) (@, <0 < @,), 


0 

so that b(O)e 0 dé ‘ K e~-O O(et), 
. 1 
6G; 


b 
ind since 4(@,) = 0, | £(6,)sinh 6, dy = 0. 
a 
Thus the contribution to the integral over Sj from 


- 


7” 





| P4 $) B sinh 6 dédy, 
oo 
taken over this range, is 
B ( $(8,)sinh 6, dy+ Ket. (6) 


In the range of integration (iii 


(6) O(1) (0, < @ < 63), 
Q 
so that | d(O)e-9 dé| < K\e-? > - O(e!) 
A; 


ind the contribution to the integral over S; from 


j 4 Ch 
| | 4 $)Bsinh 6 d0dx, 
ct 


taken over the range (iii), can be written 
B | [¢(@5)sinh 6,—¢(@,)sinh 6,|dx+Ket. (7) 


b 
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The first term in the integrand is of order «~!, but we shall see that it is 
exactly balanced by a corresponding term in the integral over S}. 


The range of integration (iv) corresponds to the surface S;. Here 


$(8) = O(1) (08, < 8 < A) 


so that 
GO, 
b(O)e-8 db| < K\e-9 2 K'\e-9 4 Ole) (A cosh—!(a/e)). 
03 
Also, since $(4,) 0 | d(0,)sinh 0, dx 0. 
i 


Thus the contribution to the integral over Sy; from 


"Lad, - 
| 54 ¢)Bsinhé@ dédx 
ot 


is B | (@,)sinh GO, dx Ke. (S) 


h 
The first term in the integrand is of order e~!, but we shall see that it is 


exactly balanced by a corresponding term in the integral over Sz. 


5. Integrals over S; and S, 


On Sz. 1, m,n 0,0, —1, and on Sz, 1, m,n 0,0, +1, so that the 


corresponding integrals are 


1 (cd | . 
| | —(h)r_ 19 : | dédn, 
| S\ cc t 0 OC\S | 
*f (1fed ofl : 
| | | —(¢)r_ 9 | dédn, 
JJ (sled/r--o 76\5 


taken over Sz. In order to examine the convergence of these integrals we 


taken over Sz, and 


< 


express 8, €, 7 in terms of z, 6, y. On the plane € = 0 


ssinh @ cos x; 


so that 
' E—2x bz coth @see x, 
n—F ztan x, 
and didn Bz* cosech?6 sec*y d6dy. 


The range of integration in the variables @ and y is 6 = 0, to 0 = 6, and 


x = bto) a, where @;(y) is the value of @ on the leading edge. i.e. 


6. eoth LS(x x)COS x Bz. 
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Over S343, 


a 6; 





"¢ 1(éd > r # * (0d | 
dEd1 Bz | sec*y d , cosech 6 dé 
sa7 X &X aes | 
Jw P\GS/ 0 ° d C¢ C=+0 
Os 
a Js 
lads _# , 4 
max . Bz sec“y dy | cosech 6 dé 
C$/(=+0 A » 
. b 6s 


nee - 19 9s 
K tan x|, —log coth $6 8, 


O(1) (2 €). 
since 
eoth AY coth| $ coth L(x x)cOos x Bz} | CAL) (¢ €), 
coth 34, coth| } sinh L| zsec x/€)| O(1) (2 €), 
und tan x , O(1) (2 €). 
- i ( | " 
\lso | (hd ‘ | dédy, taken over SJ, Sz, can be written 
wn P ( {1 * ey C | 
(p— 3) r- .9 | } dédy t- | | (d5)z 0 dédn, 
CC ~ CG Ss 


taken over the same surfaces, where ¢, denotes the value of ¢ on the 


9 


intersection of the sheet of the hyperboloid s? = e? (€—ax < 0) and the 


plane 0. i.e. d(@.) in geodesic coordinates. Now ¢d—d, vanishes when 
3 


B-(€—2) ((n—y)*+27}! (C 0) 


nd is continuous in €. if d is assumed to be continuous along the wing, 


so that in the neighbourhood of that curve 


9 


d—d, = Of B-\(E—2)4{(n—y)? +274] = O(se 


g 


). 
Therefore, applying Cauchy’s test for convergence with respect to € to the 


louble integrals 


pus | Pat | 
I3: 13, | (b—dz re | ) ag, taken over Sz, Sz, 
COU\S 

we have 
a F a = 
re a ; 

Ij .—I3. Kz | secy dy | | e? ie| — <. | sec yx ax| | e-% ie| 
i, ia € i, 6 € 
Os (« é 


ind similarly in the second case. Consec uently 
1 : 


[ [ r ds)¢ 0 (*) dédy 


taken over Sa, Sz, converge as «> 0. 
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Over Sz 


S3 
9 /] 
= [ [ ($3)¢ ox(3) dédn 
> % 
+B | dy | ($(4s)) 
b As 


_o8inh 6 dé 


t 


£B | (6(0,))¢- .9{sinh 0,-+-e~-% —cosh 0,] dx 
b , 
+B [ ($(43) )c- .9| sinh 4, —cosh 6,;| dy+ Ke, 


i 
+ B | ($(63))r-.9coshO;dxy = O(1) (z > e), 


a 


and FB [ ((43)) 


b 


c¢-.9Sinh 6, dx 

are just the terms required to balance the infinities which have already 
arisen. Thus the contributions to the integrals in (1) and (2) from the 
surface integrals over Sj , Sz, are well behaved as e tends to zero and can 
be written 


an (4) tal s)¢ 


\s\ az 


| dédn 


a 
- | (P(93) )e .o Sinh 6,—cosh 0;|dy+ Ke, (9) 
b 
the double integrals being evaluated over Sj, Sz. 


6. It now follows from equations (1) to (7) and (9), on letting « tend to 
zero, that 


P(x, Y,2) = — —_ | | rea + _ (¢—¢s), | dé&dy 4 
2 aa 
4 = | (P(95) )e .9cosh@, dx, (10) 


the double integral being taken over Sj (€ = 0); and from equations 
(1) to (4), (8), and (9), again letting e« tend to zero, that 


1/ed 2 
eo... - | { He ‘, nies <3 (hd a b3)¢ 4 dédn =. 
= 0 ¢ 


a 


+3, | ($(63))¢--p cosh 0; dx = 0, (11) 
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the double integral being taken over Sy (e€ = 0). Replacing z by —z in 
(11) and noting that the direction of integration in the line integral is 
thereby reversed, we have 


1 wae od . 
_- 2InrB | | ;( ‘I ‘i “(9 $3)¢ | dédn— 


al 
a | ($(93))¢-.9cosh0,dy = 0, (12) 
"oe 


the double integral being taken over Sj. Combining (10) and (12), we get 
two further expressions for (x, y, z): 


1 ff (ad 
X,Y, 2 ae pes léd 13 
d(x, y, 2) =A fa. Ps ” (13) 


and 


h(x, y. z) 


e , a 
zB 3? bs)y 0 dfdy + — | (P(95) )z .9 cosh 6; dx, 
b (14) 


the double integrals being again evaluated over S}. 


7. Comparison with formulae obtained by Hadamard’s method 
The formulae (13) and (14) are equivalent to those deduced by Ward (2) 
from Hadamard’s general solution. His second formula is, in our notation 
(double integrals being evaluated over S,, and an asterisk denoting Hada- 
mard’s ‘finite part’), 
1 re 


aB 








(b) 055 dédn. 


This can be written 
] a 


| © ae - 
- b—da)y — dédy — 
7B | | \P ou: +0 53 “7 7B 4 





| ((9s) )e sit 5 dédy 


93 2 if , 
—=5 | (b—$ds)¢ 05 agdn — | (0,)[sinb §,—cosh 6, | dx 
b 


a || (4—¢;)¢ . agdn+— | $(0,)oosh@, dx 
; b 


since 
a a 
in ($(45))z- .98inh 6, dx = | .98ee x dy = Kze= 
a mi 
and therefore does not contribute to the ‘finite part’ of the double integral, 
which omits terms containing negative powers of e. 
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8. Conclusion 


It is important to note that none of the formulae (10), (13), or (14), 
whether derived by Hadamard’s finite part theory, by the analogue of 
Green’s method, or by any other means, can be shown analytically to give 
the correct boundary values on the wing itself; but if we are assured of the 
existence of a supersonic potential for the problem, as we are from physical 
considerations, it must be given at points not on the wing by these 
formulae, any two of which may be regarded as an expression of the 
conditions for the existence of an analytic potential satisfying the given 
boundary conditions. Hadamard (1) states that, in general, the success 
of a direct verification of formulae obtained by Green’s type of argument 
will depend on the vanishing of the area of intersection of the characteristic 
conoid and the surface bearing the data of the problem as the vertex of 
this conoid approaches the surface; Green’s method as applied to the duly 
and non-duly inclined plane in the case of the three-dimensional hyperbolic 
equation indicates the same fact. 

Hadamard’s method has the advantage that it enables one to differen- 
tiate partially under the integral sign (even when the limits depend on 
the variables) and thus verify directly that the differential equation is 
satisfied by the formula for d but, on the other hand, such verification is 
unnecessary for the Green solution. Green’s method has the advantage 
that it is a much shorter and more direct process than the setting up and 
application of Hadamard’s analysis, which was devised for the general 
n-dimensional problem. In view of these facts, and because the three- 
dimensional problem is the one of most practical interest, it has seemed 
worth while to write the paper. The boundary value problem for oscilla- 
tory supersonic disturbance potentials can be treated in a similar manner. 
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NOTE ON A GENERALIZATION OF RAYLEIGH’S 
PRINCIPLE 
By W. J. DUNCAN (Department of Aeronautics and Fluid Mechanics, 
The University, Glasgow) 
Received 13 February 1951] 


SUMMARY 


Anv root A of the characteristic equation of a set of n linear ordinary differential 
1 i 


equations with constant coefficients and of order m is also a root of an equation of 


degree m whose coefficients are quadratic forms in the modal coordinates correspond- 


ing to A. It is shown that, when the matrices of the coefficients of each order in the 


differential equations are symmetric, the root A of the equation of mth degree is 


stationary for small deviations of the modal coordinates from their true ratios for 


mode considered. 


1. We consider a set of » linear and homogeneous ordinary differential 
equations with constant coefficients and of order m. In matrix notation 


the whole set can be written 
(A,, D™+A,,-,D"1!+...+A )r = 0, (1) 


where A, is a constant square matrix of order n, 2 is the column whose 
elements are the n unknowns, and D = d/dt. If we now assume that x is 
proportional to exp(At), equation (1) becomes 

(A,,A"+A,,_,A™-1+...+Ap)a = 0 
which may be written concisely 

A(A)x 0. (2) 

This represents n scalar equations whose eliminant is an algebraic equation 
for A, in general of degree mn. Let A, be a root (real or complex) of this 
equation and x, the corresponding modal column. Then 

(A,, Am +A,,-, A” 1+...+A)z, = 0. (3) 


In the further treatment we postulate that the matrices A,, are all sym- 
metric. When we premultiply (3) by 2, the transpose of x,, we obtain 
the following scalar equation of the mth degree for A, 


A™ x, A, 2. +A™-1x, A aXe toot i A, Z. 0, (4) 


where each coefficient is a quadratic form in the elements of x,. If 2, is 
known exactly, one of the roots of (4) will be exactly equal to A,. Suppose, 
however, that the elements of x, are in error by small quantities of the 
first order. Then it will be shown that A,, as given by equation (4), is in 


error by a small quantity of the second order. In other words, A, as given 


Ss 


[Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 1 (1952)] 
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DUNCAN 
by equation (4) is stationary for small deviations of the elements of x, from 
their true ratios. 

The above mentioned stationary property can be made the basis for 
deriving a characteristic root and the corresponding mode by successive 
approximation. Thus, from the approximate value of a root we derive an 
approximate mode and from this calculate the coefficients in equation (4). 
An improved approximation to the root can now easily be obtained from 
(4) since a first approximation is already known. A more accurate mode 
is now obtained by substituting the improved root in (2) and the whole 
process can then be repeated as often as desired. 

In the proof we require an expression in matrix form for the whole set 
of first-order partial differential coefficients of a quadratic form. Let 





Ale) = 2' Az, (5) 
where A is symmetric. Then 
OOP? x the (6) 
ex 


where the expression on the left of the equation concisely represents the 
column of partial derivatives. Hence the result of differentiating (4) par- 
tially with respect to the elements of x, can be written 


Cc A. ’ ’ ‘ - 
oe (mA 1a, Ay, Let... +2, AyxX,)+2(AMA,,+...+Ag)x, = 0, (7) 


where the multiplier of @A,/éx is a scalar. By equation (2) the second term 
in (7) vanishes and we obtain 
or, 


Ox 


0. (8) 


Rayleigh’s own theorem (1) is concerned with the degenerate second- 


order system (A, D?+A,)a = 0 (9) 


which leads to a linear equation in \?, namely 
9 ms | af 2 
Wx'A,x+2'Ayx = 0. (10) 


From the present point of view, this is not distinct from the first-order 
system 
: (A, D+A,)x = 0. (11) 
The extension to the complete second-order system was given by Duncan 
and Collar (2). 

Equation (4) has m roots and we have shown that one of these is A,, the 
root of the eliminant of (2) corresponding to the modal column z,. In 
general the remaining (m—1) roots of (4) are irrelevant and are, moreover, 
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not stationary for small variations of x,. Let » be one of these irrelevant 
roots. Then 


A(pu)a, # 0, (12) 
and equation (7) shows that 
FO. (13) 
Cx 


The irrelevant roots could be identified by use of (12). 

Suppose that A, happens to be a double root of (4). The multiplier of 
\,/éx in (7) then vanishes and we can no longer infer that A, is stationary. 
The following investigation shows that, in fact, it is not now stationary. 
Let A, be a true characteristic root and x, the corresponding modal column. 


[Then we find that (x’ +e')A(A,)(a,+e) <A(A.)e, (14) 


where € is an arbitrary column, for 
e A(A,)a, = 2, A(A,)e 


by transposition, since A(A,) is symmetric, and therefore both these 


expressions vanish by (2). Put 


€ nk (15) 

where k is a constant column and 7 is a scalar. Then (14) becomes 
(a, +k’ )A(A,)(2,+ nk) nk’ A(A,)k. (16) 
Let A= A,+6, (17) 


where we regard A, as constant and write 

flén) = (anh )AA,+E Va, + nk). (18) 
On account of (2) f(0,0) vanishes and we have by Taylor's theorem, 
when € and 7 are small, 

T(E, 0) pé+qn+ 4tré*+8En-+ biy?, (19) 
where p, g, r, 8, and ¢ are partial differential coefficients of f. This yields 
f(9, 0) = n+ 3tn’, 
und by comparison with (16) we see that g must be zero. Hence the 

equation for € when 7 is small is 


ré2+- sn +4tn? = 0. (20) 


ur 


P 
Provided that p does not vanish we must have &é of the same order as 7? so 
that t 2 
¢ a 2 
$ ai a? (21) 
<p 
in accordance with the stationary property already enunciated. But when 
pis zero, as when A, is a double root of (4), € is of the same order as » and 
the stationary property fails. The argument can obviously be extended 
to roots of higher multiplicity. 
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SUMMARY 


It is shown that, when the mass of an elastic body of linear form is supposed to 
condensed into rigid masses or particles situated at the mid-points of equal 
segments while the elastic properties of the body are unaltered, the error in the 
frequency of any natural mode of oscillation ultimately varies inversely as the 
square of the number of segments. When the masses are not placed at the mid- 
yoints of the segments the error varies ultimately as the inverse first power of the 
imber of segments. The general conclusions are checked by numerical examples 
ind by comparison with cases where exact solutions for the segmented bodies are 


\ ailabl . 


1. Introduction and conclusions 

[x a widely used method for the approximate treatment of the dynamics 
of bodies which are both massive and elastic the mass is supposed condensed 
into a finite number of particles (or sometimes extended rigid bodies) 
embedded in an ideal massless substance possessing the same elastic proper- 
ties as the body simulated. In this way a body having infinitely many 
degrees of freedom is represented by a dynamical system having a finite 
number of degrees of freedom. Hitherto this method has been applied 
somewhat uncritically, and the aim of the present paper is to examine the 
errors inherent in the method and to seek the means to minimize them. 

My interest in this question was aroused when I observed that the error 
in the fundamental frequency of torsional oscillation of a uniform cantilever 
appeared to vary inversely as the square of the number of equal segments 
into which it was imagined to be divided and there was an indication that 
the same rule applied also to the first overtone frequency (1). A proof that 
this rule is indeed correct in the limit when the number of segments tends to 
infinity is supplied in the present paper; it is shown that the rule holds for 
ul the natural torsional frequencies of the cantilever. The proof depends 
on an exact solution of the dynamical equations of the segmented body and 
it can be extended to other conditions of support and to other physical 
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problems, such as the longitudinal oscillations of uniform rods and the 
transverse oscillations of uniform cords in tension. In all these cases the 
particles are placed at the middles of the segments, which are of equal length. 

[It is also shown that the inverse square law of error holds in the limit for 
any uniform or non-uniform ‘linear’ body, provided that the particles are 
placed at the mid-points cf segments of equal length. This can be deduced 
from Rayleigh’s Principle, provided that the errors in the modal ratios 
measured at the particles are either zero or tend to zero in proportion to 1/n 
or some higher inverse power of n in the limit. An argument is advanced 
which leads to the conclusion that the modal errors in fact tend to zero in 
proportion to the inverse square of the number of segments. The generality 
of the inverse square law of error in the frequency is thus established for 
‘linear’ bodies. As a numerical example, the errors in the frequencies for the 
fundamental and first two overtone modes of a uniform cantilever beam 
performing flexural oscillations are obtained for numbers of segments up 
to 10 and the correctness of the inverse square law is verified. No exact 
solution for this problem with an arbitrary number of segments is known. 

In all the cases investigated the proportional error in the frequency, with 
a given number of segments, is substantially higher for the first overtone 
than for the fundamental and, generally, increases rapidly as we move up 
the sequence of overtones. The rough rule that the use of 13 particles per 
wave-length ensures that the error in the frequency is not greater than 
about 1 per cent. fits the results hitherto obtained. 

[f the particles are not placed at the mid-points of the segments the errors 
are greatly increased and the error in the frequency is ultimately propor- 
tional to the inverse first power of the number of segments. 

Knowledge of the manner in which the error in the frequency depends 
on the number of segments can be used to improve the approximation when 
results have been obtained for two or more numbers of segments, in accor- 
dance with the procedure of ‘the deferred approach to the limit’ (2, 3, 4). 
This is illustrated by a numerical example. 


Acknowledgement. I am indebted to Mrs. 8. D. Silvey for assistance with 
the calculations discussed in section 4 of this paper. 


2. The error law for some exactly soluble problems 


In some simple instances it is possible to obtain an exact solution of the 
equations of free oscillation of a body having its mass condensed into equal 
particles or rigid elements situated at the mid-points of n equal segments 
when » is arbitrary. Since we are not primarily concerned with obtaining 


these solutions, and since they can easily be verified, we content ourselves 
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he with quoting the results and remarking that they are most readily estab- 
he lished by finite difference methods.t 

h. The following three problems have an identical mathematical formula- 

‘or tion, with suitable interpretation of the symbols: 

re 1) The small transverse oscillation of a uniform and perfectly flexible 
ed cord in tension. 

0S 2) The longitudinal oscillation of a uniform and straight thin elastic rod. 
n 3) The torsional oscillation of a uniform elastic shaft. 


ed | For definiteness we shall discuss problem (3), but the scheme of re-inter- 


in pretation of the symbols for the other problems is given in Table 1. The 
ty ndition at either end of the body will be end fixed (@, uw, or y zero) or 
for end free (zero applied couple, longitudinal force, or lateral force). 
he In accordance with what has been said, we take a uniform elastic and 
mM wsive shaft and suppose it divided into n segments of equal length. A 
up | flywheel whose moment of inertia is equal to that of the whole segment is 
uct laced at its mid-point and the elastic reaction to its displacement is the 
mn same as for the actual shaft at the same position. First we take the case 
ith | of the cantilever shaft (one end fixed and the other free). The fixed end is 
me the origin from which the abscissa 2 is measured. The exact value of the 
up non-dimensional frequency parameter as defined in Table 1 for the rth mode 
e] f free oscillation or (r—1)th overtone is 
an 

; 2n2| ada 2r—1)z (2.1) 

l 2n } 

OTS i 
or vhile the angular amplitude at the sth flywheel in this mode is 


| 2s—1)(2 l)a 
_= asin 1 ld ad (2.2) 
ids 


tn 


en , ; ; ' . . 
| where ais an arbitrary constant. Now the displacement in the corresponding 
free oscillation of the continuous shaft is 
| 


¢ 6sin : (2.3) 


ith , ; ; 
vhere 6 is an arbitrary constant. At the sth flywheel the abscissa is 
2s—1)l/(2n), so 
~ . (2s—1)(2r—1)z 
G b sin : (2.4) 
47 
the nd this agrees exactly with (2.2) if we make b = a. Consequently there is 


ual | exact agreement of the modes of displacement at the flywheels. Between 


nts . , — . . 
See, for example, Karman and Biot, Mathematical Methods in Engineering, chap. xi, 
ing section 6. The problem of the oscillatio f a massless cord loaded with equal and equally 
ed particles was completely solved by Lagrange, who claimed no more than to have 
ves 


tematized the work of his predece Ssors 
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TABLE | 


Corresponding Physical Quantities, with Symbols 


TRANSVERSE 


TORSIONAL LONGITUDINAL OSCILLATION OF 
OSCILLATION OF OSCILLATION OF UNIFORM CORD 
UNIFORM SHAFT UNIFORM ROD IN TENSION 





Anqular displacement Longitudinal displacement | Lateral displacement 
g I { 7 I 


6 u Y 
Corresponding Moment of inertia Mass of unit length Mass of unit length 
basi quantitie s of unit le ngth m m 
J 
Torsional stiffness Force per unit Tension in cord 
(Moment per radian) extensional strain = 
of unit length AE 
( with 
A Area of section 
E Young's modulus 
Non-dimensional Jl?" ml*w* ml? a" 
X x x 
frequency parameter C AE i 
l length of body. w 27 times frequency. 


the flywheels the displacement varies linearly with the abscissa and the 
displacements of the segmented and continuous bodies do not agree. 

For a given mode r the argument of the cosine in (2.1) tends to zero as ” 
tends to infinity. If we retain the first three terms in the expansion of the 


cosine we obtain the approximation 


(2r—1)®x2 (2r—1)4xr4 


Qn 
Xnn —_ (2.5) 
i 4 192n2 
: (27 — 1)?77? . 
and lim «,. (2.6) 


rn 
no 4 


which is the true value of «, for the continuous shaft. Hence we have the 


error 
(27—1)47r4 r ; ; — 
Ern X-— Nyy -_—-| higher inverse powers of 7. (2.7) 
192n7 
and in the limit . " ie 
Nn, (2r— 1)*2r* . 
(2.3) 
¥ 48 


This shows that the error in the non-dimensional frequency parameter for 


any given mode ultimately varies inversely as the square of the number of 


segments and that the proportional error for a given n increases rapidly as 
we move up the scale of overtones. The convergence of en? for the funda- 
mental and first overtone is shown in Tables 2 and 3 respectively. Especially 
for the fundamental, the convergence is rapid. 
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Let us next consider the shaft with both ends fixed. Here we find that the 


exact value of the frequency parameter of the segmented shaft for the rth 


mode is 
Yr 
\ 2n-| 1— cos —}, (2.9) 
Nn 


while the displacement at the sth flywheel is 


r(2s ] )ar 


2n 


A asin (2.10) 
und there is again exact agreement of the modes of the continuous and 

segmented shafts at the flywheels. From (2.9) we deduce in the same manner 
g before that 


(2.11) 


vhere ¥ YT". (2.12) 


‘he general conclusions are the same as for the cantilever shaft and numeri- 
ul calculations for the fundamental and first overtone show rapid conver 


sence of n7e to its limiting value. 


For the shaft with both ends free we find that a,,, is again given by (2.9) and 
by (2.12), so (2.11) is again valid. The modes of oscillation are given by 
r(2s—1)z 9 
0 a cos (2.13) 

2n 
ind there is again exact agreement of the modes for the continuous and 


St smented shafts at the fly wheels. 


3. General consideration of the law of error 
For definiteness we shall consider the torsional oscillation of a shaft. in 
general non-uniform: we shall suppose the shaft to be divided into n seg- 


ents of equal length and the whole moment of inertia of each segment will 


TABLE 2 
Convergence of en® for the Fundamental Mode of a Uniform 
y atale ver Shaft 
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TABLE 3 
Convergence of en? for the First Overtone Mode of a Uniform 
Cantilever Shaft 


Yr 
€ x Fr 
Y ! en” 
I 
13°6569 8-5497 34°20 
tio 2066 37°86 
| 19°7541 2°4525 39°24 
5 20°6107 1°5959 39°90 
I ee ee 41°09 


be supposed condensed into a rigid flywheel situated at the mid-point of 


the segment. The elastic couples are the same as for the actual shaft in the 
same circumstances. Other cases of the oscillations of bodies of linear form 
can be treated by methods similar to those used here. 

We shall base the argument on an assumption about the natural modes 
of oscillation of the segmented shaft which we shall later examine and 
justify 

The departure of the modal ratios of the segmented shaft, as measured 
at the flywheels, from those for the corresponding mode of the actual 
shaft are either zero or tend to zero as ” tends to infinity ultimately in pro- 
portion to the inverse first power of x or some higher inverse power of n. 
The possibility that the modal errors may be zero has been demonstrated 
in section 2. 

We shall calculate the frequency of oscillation of the segmented shaft 
in any pure mode by equating the kinetic energy at mid-swing to the 
potential energy at the end of a swing. Thus, if the potential energy is V 
and the kinetic energy is w?7', where T and V can be found when the mode 


of oscillation is known. we have 
w : (3.1) 


In accordance with Rayleigh’s Principle, if we apply this formula with 
7’ and V calculated for modal ratios which are in error by small quantities 
of the first order, the error in w is of the second order. Let us then adjust 
the mode for the segmented shaft so that it agrees exactly with that for the 
actual shaft at the flywheels. In accordance with our assumption about the 
modes, this either leaves the value of w for the segmented shaft unaltered 
or changes it by an amount ultimately proportional to n-? or some higher 
inverse power of n. 

We can now pass on to calculate the difference in the values of w? for the 


actual shaft and for the segmented shaft with its adjusted modal ratios. 
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First let us consider the difference in the values of 7' for the actual and 
segmented shafts when the displacements agree at the flywheels. Let us 
take a segment of length 2h whose mid-point is at « = €. Since the length 


f the shaft is / and there are n segments we have 


l 
h ; (3.2) 
2n 
Let the polar moment of inertia per unit length now be J(x), let the mode 
of displacement of the actual shaft in the mode considered be 
0 = f(a), (3.3) 
und let F(a Sar. (3.4) 


Then the contribution of the seement to 7' for the actual shaft is 
T, =4 | J(x)F(a) da. (3. 


—h 


oso 
or 


We shall expand J(x) F(x) by Taylor’s theorem as far as the second power 


of (x—€&) and so evaluate the 


ntegral approximately. In this way we get 


T, = hd (E)F(E)+- RAST (€) P(E) +25 (E) F'(€E)}+- J" (E) F(O)}, (3.6) 


t 


where the error is of the order h°. Next, the moment of inertia of the flywheel 


epresenting the se oment is 
R | J(x) da 2h J(€) th? J”’(E) (3.7) 


to the same order of approximation as before. Thus the contribution of the 


flywheel to 7’ for the segmented shaft is 
‘ly LI F( hJ(€)F(E)+ Fd" (€) F(E). (3.8) 

Hence the difference for the segment is 
T_—T. hs S(E)F'(E)+-2S'(E) F'(E)}. (5.9) 


By addition of the results for all the n segments and use of (3.2) we see that 
the error in 7' for the segmented shaft varies ultimately as n-*. It is worthy 
of remark that if we had taken the moment of inertia of the flywheel to be 
2hJ(€). i.e. to correspond to that of a uniform cylinder whose length is that 
if the segment and whose cross-section is that at the middle of the segment. 
we should still find the error varying ultimately as n-*. 

Next we have to consider the error in V for the segmented shaft. Between 
an adjacent pair of flywheels the shaft is unloaded and its energy is the mini 


mum consistent with the displacements of the flywheels. The displacement 
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of the actual shaft is the same as that of the flywheels at the sites of the fly- 
wheels and its additional displacement between the flywheels is the same as 
for the segment of shaft with fixed ends and subject to the actual inertia 
loads. Further, the addition to the potential energy of the segment is just 
that of this displacement of the segment with fixed ends. Now as we vary h, 
with / small, the inertial load is proportional to h and the flexibility of the 
segment (regarded as having fixed ends) is also proportional to h. Hence the 
relative displacement is proportional to h? and the elastic energy (which is 
proportional to the product of the load and the displacement) is proportional 
to h’. Thus the defect of elastic energy of the segmented shaft per segment 
is proportional to h3 and the total defect ultimately varies as n~?. 

Since the defects of 7’ and of V for the segmented shaft are both ultimately 
proportional to n~*, it follows from (3.1) that the error in w? ultimately 
varies in the same way. Since the error in w? introduced by adjusting the 
mode of the segmented shaft varies ultimately as n-? or some higher inverse 
power of n, we conclude that the error in the frequency of any given free 
mode as given by the segmented shaft ultimately varies inversely as the 
square of the number of segments. 

In order to complete this investigation we have to justify our initial 
assumption about the errors in the modal ratios at the flywheels. To do this 
we shall adopt the ‘Lagrangian’ method and the same finite number N of 
displacement functions and corresponding generalized coordinates for the 
actual and segmented bodies; finally we shall make N tend to infinity. The 
generalized coordinates may conveniently be chosen as the normal coordi- 
nates of the actual body. 

Since the elastic specification of the segmented body is the same as that 
of the actual body, the two sets of stiffness coefticients in the dynamical 
equations are identical. The typical inertia coefficient A,, for the actual 
body is given by 1 
A, | J (x) f(x) f(a) da, (3.10) 

0 
where /,(v) and f,(v) are the displacement functions corresponding to the 
generalized coordinates q, and q, respectively. By a treatment altogether 
similar to that which leads to equation (3.9) we deduce that the difference 
between A,, and the corresponding coefficient for the segmented body varies 


rs 


9 


ultimately in proportion to n-*, when n tends to infinity. We then deduce 
that the difference in the modes must vary with » in the same way and our 
assumption is justified. Since the modal adjustment is of the order n~*, the 
change in frequency consequent on this alone is of order n-* and is thus 


negligible. We could deduce immediately from the dynamical equations 


that the frequency difference is of order n-?. 
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ly. | 4. An illustrative numerical example 
as | As an illustrative example we take the flexural oscillations of a thin 
tia | uniform cantilever beam and make the usual assumptions of the Bernoulli 
ist | Ruler theory. The non-dimensional frequency parameter « is given by 
} 
1.2 
ml W 
he | (4.1) 
El 
he 
st where m,/, and w have the meanings given in Table | while /'/ is the constant 
‘al flexural rigidity of the beam. The true values of « are the fourth powers of 
nt the roots of the equation 
cos cosh x i} 0. (4.2) 
‘ly [here is no known exact solution for the oscillations of the segmented beam 
‘ly when the number x of equal segments is general. The frequencies and modes 
he rv each value of x were obtained by the matrix iteration method (5). The 
rs iss of each segment is supposed to be condensed into a particle situated 
ree tits mid-point. It will be seen from Table 4 that the calculations confirm 
he that the errors in the non-dimensional frequency parameter for the modes 
f oscillation examined accord with the inverse square law. The results for 
ial he errors in the modes are also in conformity with this law. 
his a ; _— 
5. Effect of placing the masses away from the centres of the 
oO - 
) segments 
1e ° ° . ° . . . 
: It is evident from the investigation in section 4 that we should find a term 
he -" - mea 
, proportional to h? in (7',—T7.) if the flywheel were not placed at the centre 
Cl : . . mL 
f the segment and we should consequently find the error in 7’ for the whole 
shaft tending to zero in proportion to n-! in the limit. Consequently we 
lat ; a 
) should expect the error e in the frequency parameter to tend to zero in 
Ca rms: > 
| cordance with the law that en tends towards a constant. This is confirmed 
la 
TABLE 4 
10) Flexural Oscillations of Uniform Cantilever 
Dependence of Error in Frequency Parameter on Number of Segments 
he ml4u)? 
1€1 El 
1ce Sé ; ton odeé 
1es8 x 455°522 l é€ value x 3500°55 
ice € é é en? n r en® 
yur ; 
he 55 { 2 1,Q12°05 : 

2 2,229°95 20,069°55 
tus } 588-77 9,420°32 
ns 22 17 507°70 5 390°990 9,922°50 

9,950°0 
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by the results for some cases where an exact solution for the segmented body 
(with a general number of segments) can be found. 

Let us take the uniform cantilever shaft considered in section 2 and now 
place the flywheels at the outer ends of the segments. It can be shown that 
the exact value of the non-dimensional frequency parameter in the rth mode 
for the shaft with n segments is 
‘ (27—1)a e 
oe 2n*| 1—cos ; (5.1) 

2n+1) 
When we retain, as before, the first three terms in the expansion of the 


cosine we get . eo e 
n*(2r—1)?x7?2  _n?(27— 1) 4214 


77 


(2n+1)? 12(2n+-1)4 * 


Evidently lim «, 
which is the true value of «, for the continuous shaft. But we find now 


E.,, X,—,, —> (5.4) 


when v is large. The convergence of ne for the fundamental mode is shown 
in Table 5. 


TABLE 5 
Representation of Uniform Cantilever Shaft by Segmented Body 
(masses at outer ends of segments) 
Convergence of ne for Fundamental Mod: 


€ He 
1°467 1°467 
0°939 I°d79 
0°0355 54 
} 0°537 150 
0°442 ite) 
O°2334 4 


In order to obtain the fundamental frequency correct to 1 per cent. we should 
have to use 50 segments. Comparison of Tables 2 and 5 shows the immense 
loss of accuracy when the masses are moved away from the centres of the 
seoments. 

[t can be shown that the amplitude at the sth flywheel in the rth mode 
is now 


s(2r—1)a 


6, asin- (5.5) 
2n+ 1) 











This 


corre 


whic 
Tl 
segn 
the : 
exce 
of t« 
CU; 
part 
indi 
6. I 
L 
the 
and 


fror 


req 


act 
len 
of 


ha 








ody 
LOW 


hat 


ode 


wn 





wuld 
nse 


the 


ode 








REPRESENTATION OF MASSIVE AND ELASTIC BODIES 107 


This flywheel is situated at the distance s//n from the root, and from (2.3) the 
corresponding displacement of the continuous shaft is 


A a nd (5.6) 
2n 
which cannot be made to agree with (5.5). 

The value of «,,, when the flywheels are placed at the inner ends of the 
segments can be obtained from (5.1) by changing (2n+-1) into (2n—1) in 
the argument of the cosine. Hence the results are qualitatively as before 
except that the frequencies for the segmented body are all too high instead 
ot too low. 

Calculations on the flexural oscillations of a uniform cantilever with the 
particles placed at the outer or inner ends of the segments also yield results 
indicating that en tends towards a constant. 

6. Use of the law of error to improve the approximation 

Let a, be the value of « for the body with n segments and suppose that 

the error is proportional to the inverse pth power of n. Then if we know p 


ind have calculated «,, and «,, we may find « with improved approximation 


irom the equations ' v-t-km P| 
. (6.1) 

vt+kn-P? } 
, : , x, 4 X, Mm” pa 
which yield . (6.2) 


nr’ m” 


If we regard p also as unknown we could find it and then the value of « from 


1 set of three equations similar to (6.1). 

\s an example let us take the results from Table 2 for 4 and 5 segments 
nd assume p 2. Then equation (6.2) yields a 2-46728. The error in 
this (to 5 places of decimals) is 0-00012, whereas the error in a; is 0-02023, 


e. about 169 times as great (see also 1). 


7. Rough rule for the number of masses required 

Let the masses be placed at the middles of the segments and let it be 
required to find the frequency in some particular mode of oscillation correct 
to 1 per cent. Then the permissible error in the frequency parameter (which 
is proportional to the square of the frequency) is 2 per cent. A survey of the 
results quoted in this paper and some others suggests that this standard of 
vecuracy will be reached when the number of segments per complete wave- 
length is about 13. This refers to the frequency uncorrected by the methods 


ot section 6 
8. An electrical analogy 
It was pointed out by Heaviside that the inductance of a telephone cable 


has to satisfy a certain relationship in order that transmission of speech 
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should be ‘distortionless’ and, in practice, the satisfaction of this condition 
requires an increase of the inductance. It was at first proposed by Krarup 
to achieve this increase by lapping the cable with thin wire of high per- 
meability (continuous inductive loading), but Pupin showed that the same 
result could be more conveniently obtained by inserting inductance coils at 
regular and not too large intervals (discontinuous inductive loading). The 
rule for the spacing of the coils given by Pupin is as follows (the quotation 
is from Fleming (6)): 

‘If there be a non-uniform cable line loaded with inductance coils at equal inter- 
vals, and if we consider the total inductance and resistance to be smoothly distributed 
along the line, then these two lines, the non-uniform and uniform lines, having the 
same total resistance and inductance, will be electrically equivalent for transmission 
purposes as long as one half of the distance between two adjacent coils expressed as 
a fraction of 27 taken as the wave length, is an angle so small that its sine has 
practically the same numerical value as that angle in circular measure.’ 

This rule is in substantial agreement with that enunciated in section 7 
above. 


9. Questions to be examined later 
[t is hoped to extend the investigation of the segmental representation 
of bodies to those having two- and three-dimensional extension and to cover 


their mechanical admittances as well as their natural frequencies. 
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{ METHOD FOR THE EVALUATION 


er- 
me . . n 
xm (2 7) | exp( se) at dx 
a 3 
' i) 
he 
on a Sie . ; 
By H. J. GODWIN (University College, Swansea) 
er ~ ~17 
Received 5 June 1951] 
¢ 
he 
ion SUMMARY 
AS ‘ 
las The integral | | (2/77) | exp $12) dt) dx is evaluated by expressing it as a series 
0 ri 
of reciprocals of factorial functions of n; the accuracy of this method increases with n- 
t is also shown how the method may be modified to improve the accuracy for small 
s of 
on 1. Introduction 
er L 
1 We denote \ (2) z)exp(— 44?) dt by g(x), and | x™{g(x)}" dx by I,,,. These 
0 
integrals arise in certain problems connected with the normal curve of 
'o error; for example, /, , can be used to obtain mean values of rank differences 
in random samples, and J, ,, (” 0) to obtain variances of such differences 
eS see Godwin, 1). The purpose of this paper is to express J, ,, a8 a series of 
wa reciprocals of factorial functions of n, and to discuss the accuracy of the 
] evaluation of the integrals in this way. 


2. Derivation of series for / 


Define a sequence of polynomials in x by the equations 


d 
} P(m.r.2x) rzP(m,r—1,2)-4 P(m,r—1,2z) (1) 
ax 
and P(m 0.7) 7 (2) 
Then, from (1), we have 
_— d oii ‘ 
exp(Sra*)P(m,r, x) ; fexp(Sr2?)P(m,r i205. (3) 
ar 
da ; 
Also (47)! exp($2?). (4) 
dg(x) = i 


[Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 1 (1952)] 
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Now 
0 
ie P(m, 0, x)g" (477) exp(42) dg 
i 
q'' 1 g=0 
(377) exp(da®) Pom, 0,2)| 
n+] —_ 
0 
titi . 
| 7 (a7)exp(2*) P(m, lx) dg, 
i 
1 


on integrating by parts and using (3). Integrating by parts k times we have 


k-1 ' 
" Mt. 
T= S (An) P(m, 7, 0) R(m, k, n), (5) 
roart rm (n+r-+1)! 
7 { 
where 
0 
fe ! 
n! 
Rim. k.n) (Ar) ht Dgn th -exp{ i(k L)v*! P(m. ki) dg. (6) 
al = (7 k)! ~ 
1 


Since all terms in (5), including R(m.k,n), are positive, R(m,k,n) > 0 
as k-> o. Also, from (6), 
n+] 
Rim, k,n-+-1) Rim. k,n). (7) 
nt+kh+] 
The series in (5) gives the required evaluation of J,, ,,; the coefficients may 
be calculated by first obtaining the polynomials P(m,r,a) from (1) and (2), 
but an alternative method, and one which gives information about the 
asymptotic behaviour of the coefficients, is as follows. 
Define a function ¢,, by 


2™ exp(4z?) bn | exp( 12) ai). (8) 


0 


Then ¢,,(w) is an even or odd function of w according as m is even or odd, 


and has a singularity at w = (47)! (and thus also at w (47)'). Also, 


> 


for w near (477) . 


d,, (Ww) {(47) wt—[log{ ($a)! —w}-tyhim Db (w), (9) 
where %,,(w) is regular at w = (47)! and ¢,,f(47)!} = 26m-D, If 
Im\2 #) exp( 12°)4,,(texp 422) 4 exp( 522) de} (10) 
: 0 


then f,,(z, ¢) is a regular function of z and ¢ in some domain surrounding the 
origin, and 


xt 


Sn(2,t) = > t*y,,(z) (11) 
0 


n 
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in this domain; in particular y,(z) 2”. From (10), 


f(z t) +2f»(Z, t) (1—2t)< f(z t) 

" ot’ 

whence, substituting from (11) and equating coefficients of powers of ¢, 
ve have 


- 
Xn\*) 2Xn—1\) Xn 1(2). 


ymparison of this with (1) gives 


P(m,r.z) 
XZ ‘ (12) 
f} 
On putting z 0 and ¢t win (10), and using (11) and (12), we have 
~ P(m,r,0)w" 
d (7) * ‘9 e (13) 
si Lod r! 


0 


z)\(c—2z)-Hlog(« z)-1* and f(z) » a, 2", then 


l \« - = 


Now it has been shown by Eggleston (2) that if f(z) has a singularity of the 


a, ~ d(c)(log n)Fe— +, 


Eggleston actually proves this under the assumption that / is an integer, 
ut the proof can be extended at once to all real k.) Combining the two 


singularities of d..(w) and using (9) and (13) we have 


P(m,r, 0) 
esis ~ 2 (log r)i Uf ,/(477)} (r+1) (14) 


' 
} 


miseven. (Ifr—m is odd then, from the remark following (8), P(m, 7, 0) 
8S ZETO. } 


From (13) we have 


P(m,7r,0) l [ 4,,(w) _ 
dw, (15) 
) 271 wrt 
( 
vhere Cis a contour surrounding the origin, but not enclosing w + ($2r). 
If in (15) we put wu | exp(—4C*) dZ, then, from (8), 
0 
. . (r+1) 
P(m r.U) | 172 , > 
, = 2 | exp(— $C?) dZ dz, (16) 
/ = 77 t 


| 0 
where I’ is a contour surrounding the origin in the z-plane. Consequently 
P(m,r,0) is r! times the coefficient of 2’-” in the expansion of 
= r+] 
” . 
| exp(— 4°) dé (1—$2?+- A24—...) +», 


« 
0 
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P(m,m--2yv, 0) 


If we put Yn To (dar) hont2v +9) (17) 
(m+ 2p)! 
: ' 9,,\! 
: . “ n!(m+2v)! 
then, from (5), ies » Ym.y 5 -. (18) 
pa (n+m+2v+1)! 
y=C 
From (14) and (17), 
ay — Pim DSlog(m | 2y) Len 1) (19) 
fm, = ( So = j . v) 


3. Improvement of accuracy for small » 

In this section, to fix ideas, we consider the case m = 0 only, but for other 
values of m the method is similar. The terms of (18) decrease rapidly at 
first, but the ratio of consecutive terms tends to one, and the series ulti- 
mately converges slowly. To obtain a more accurate answer one could 
substitute expansions for the y,,, and sum the resulting infinite series. The 
presence of the logarithm terms in y,,, would, however, make this im- 
practicable. (If m 1, the first term in the expansion of y,,, contains no 
logarithm, but subsequent terms do.) The integrals Jy ,. Joo, Jo3, and I), 
can be evaluated explicitly in terms of elementary functions and it is found 
that for these greater accuracy is obtainable by fitting an expression 

“y = + » 

— ats 2v(2v—1) © m) 
to the last few known y,,. The series arising can then be summed by 
algebraic methods; for 


2N+a 2N+a+1' 2N+a42 
(—1)?-* 20-1 | —_ | : | 
2N Lab l 2N 1 h 


and, on using Euler’s method for the transformation of slowly convergent 
alternating series, 


l I I = t! 
nm wi‘ a+2 x 2 In(m—+1)...(n-+#) 
0 


For /,, the remainder after the term in y, ,, (i.e. R(0, 37,4) in the notation 
of (5)) is 7955 x 10-1. The use of expression (20) with three terms reduces 
this to —1510-! and with four terms to —610-!°, The accuracy 
obtainable by this method will increase with n, since the terms with larger v, 
which are fitted less accurately by (20), become relatively less important 


as n increases. 
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The particular forms of equation (20) which were used in the preparation 


of the tables at the end of the paper, and their sources, are as follows. 


Yo16: Yo.17> aNd Yo1 Satisfy 


4-198376 28-573103 
0-604468 “ . 


Yo, a mg 
2y 2v(2v—1) 
Yo.15> Yo0.16> 0.17 and Yoig Satisty 
e _., . 0378712 72-762108 455-094355 
Vo.: 0-590526 - + A 
2v 2vn(2v—1) © 2v(2v—1)(2v—2) 
¥19 and yj, 1» Satisfy 
: 0-925764 
Via 1-803484 — 
2v+-1 
1g: Vig» and vy, 4 satisfy 


_ ss 1-478476  5-101956 
Vv 1-S17656 : 
ss 2v+1 2v(2v+-1) 
99 ANd Yo 4 Satisfy 


16-123777 


Vo, 4-716121 
" 2v+2 
Vos: Yo9 ANd yo 19 Satisfy 
. 28-337187 118-857330 
5-014009 } a 
2v+2 (2v-+-1)(2v4 2) 


4. Tables 


\s examples of the cases in which y,, , tends to zero, to a constant, or to 


st 

infinity as v tends to infinity, the values of y,, , and J,, ,, have been evaluated 
form = 0, 1, 2 and various values of v and n. The tables of coefficients y 
which are given, and the expansions at the end of section 3, enable J, ;, 
I,,, and J, to be evaluated to ten, ten, and nine places of decimals respec- 
tively; thereafter the number of places obtainable increases by about one 
for each increase of one in n. (For the smaller values of n, direct methods of 


evaluation are possible.) 


5. Extensions of the method 


f(x){g(a)\" dx may be evaluated by expanding f(x) in a power series and 


1 

integrating term by term or by defining functions from f(x) analogously to 
the way in which P(m,7r,2) is defined from x2” and so obtaining a series of 
for the 


’s similar to the y,,,. If an expression involves a number of J 


5092.17 


m,n 


I 
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same m and different n, then it may be easier to obtain the solution in the 
form > « 
pa 
y 


substitute these; this is the case, for example, in the statistical application 


before substituting the y’s, than to evaluate the J, and 


» ay 
véiliv / m,n 


mentioned in section 1. 
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TABLES 








n lo.n lien | Io.» 
I *79788 45608 *50000 00000 | *53192 3041 
2 *49735 99545 "15169 O1138 *12353 O105 ‘ 
3 *33490 29304 09939 923609 "04945 9562 
4 *26205 22503 "00035 10707 *O251I 1554 2. H. 
5 *21569 24955 "04153 505061 "01400 0372 
6 *18344 12341 | "03041 74237 ‘00928 6087 | 
7 *15907 36772 "02327 25892 | "00625 7560 
8 *I414I 01257 "01539 95154 "00446 3387 
9 *12692 61705 "01492 23808 "00328 6807 
10 *II515 25543 "01235 18558 "00249 2696 
| 
II "10539 009044 "01039 65415 "00193 0666 
12 "0Og710 15794 "00597 39045 "OO153 5301 
13 *OQOI3 12162 | "00706 43175 "00123 5262 
14 | "08405 34996 | ‘00665 73679 “OOIOI 3459 
15 07574 06723 "00555 07355 "00054 0222 
16 "07407 239058 "00522 22750 "00070 4463 
17 "06992 360222 "004606 46756 "00059 6554 
18 ‘00621 62441 "00419 21232 "00050 9652 
19 "06288 32145 "00375 S182 "00043 5945 
20 } "05997 04410 "00343 99534 "00038 0755 ; 
v 0 l Yo 
Oo 1°25331 41373 15500 10 *74312 3679 
I "99435 0060210 07051 II *73542 9033 
2 *QOIQ5 53054 05054 12 *72055 9707 ; 
3 85682 23651 79658 13 °72236 8784 
4 "82680 09517 98864 14 "71674 4288 
5 15 "71159 9074 
6 16 *70686 4098 
7 17 "70248 3802 
8 | *76185 53515 54 18 ‘69841 2868 
9 75184 12218 35 ' 
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n the on . 


and : — 
8) ; 4897 1°Q6570 12432 15302 
ation I +t pO S< 2°57702 39011 05953 
2 92503 50335 
3 3°17170 O1003 
4 3°35700 15885 
5 ] s*SOSOI 55005 
some 6 52977 ~ a ig 
. 7 3°73575 C5104 
AT101 x 1044 382814 36967 
9g } 3°90993 23143 
10 3°98322 22110 
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\ TABLE OF THE GENERALIZED RIEMANN ZETA 
FUNCTION IN A PARTICULAR CASE 


By E. O. POWELL 
(Microbiological Research Dept., Ministry of Supply, Porton, Wilts.) 
[Received 20 March 1951] 
SUMMARY 
A ten-decimal table of the generalized zeta function, f(s, a), is presented, with 


8 } and a 1-00 (0-01) 2-00 (0-02) 5-00 (0-05) 10-00. 


Modified second central differences are provided. 


Introduction 
THE occasion for computing the present table arose out of a problem 
involving a function which should satisfy, among other conditions, the 
difference equation 


f(a)—a+ = f(a+1). (I) 
Equation (I) is seen to be satisfied formally by the divergent series 
¥ (n--a) 
. 


an expression which suggests the application of the generalized Riemann 
zeta function due to Hurwitz (1). This function is defined for re(a) > 0 by 


{(s,a) = } (n+a)-* (re(s) > 1) 
0 
or by 


l V1-s 
C(s, a) om Ss — | (re(s) > 0,s 4 1), (II) 
7 ar (n +-2)* | s| 


or more generally for all except positive integer values of s by 
'(l—s) ( (—wa)*-1te-2" du 


2m . i—¢-" 


c 


C(s, a) (ITT) 

where the contour proceeds from -+00 on the real axis towards the origin, 

encircles the latter anti-clockwise, and returns to +-00, 

the poles u }+-27ni of the integrand. Hence it is easily established that 
C(s,a)—C(s,a+1) a-s, 

and so the general solution of equation ([) may be taken to be 


C(4,a) 


apart from the ‘arbitrary’ periodic term depending on initial or other 


conditions. 


(Quart. Journ. Mech. and Applied Math., Vol. V, Pt. 1 (1952)] 
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A TABLE OF THE GENERALIZED RIEMANN ZETA FUNCTION 117 
The generalized zeta function does not appear to have been tabulated 
except in the special cases where it coincides with other functions; for 


example 


d2 
and C(2, a) flog ['(a)}. 


da* 
The notation used here is that of Whittaker and Watson (2). Mitchell (3) 
discusses a still more general function which he calls (a, s|z) defined by 


r(l—s) f (—w)*-te-4* du 
f(a, s|z) 2 ; 


or with suitable restrictions on the arguments by 


x 


C(a, 8 \2 S 


hunt (n t ays’ 
0 


an 


so that 


(a, s|1) (Mitchell (3)) = ¢(s,a) (Whittaker and Watson (2)). 


Properties of the generalized zeta function 
(i) By applying the Euler—Maclaurin summation formula to the right- 
hand side of (If) an asymptotic expansion convenient for actual computa- 


tion is obtained: 


B,  B1.3.5 , Bg1.3.5.7.9 


‘(t,a)~ 2a la | i a wns y 
o\2 - 2! 9q3/2 4! 23q7/2 6! 25q11/2 (IV ) 
where the B, are Bernoulli’s numbers 1/6, 1/30, 1/42, ete. 

(ii) An ascending series in powers of a, convergent for |a| < 1, may be 


obtained from equation (II1) by writing (4, 1+-a@) as 


(4) [ | w)—te—YNe-au dy 


expanding the factor e-*” and integrating term by term: 
: « (2n)! ./2n+-1 . 
¢(4, 1+a) 7 , $1 ——]}(—@)*. (V) 


rhe coefficients here diminish only very slowly, so that the useful range 
of the series is limited. 

(iii) In common with the solutions of other difference equations of the 
general type of (I), ¢(4,@) possesses a finite multiplication theorem which 
can be obtained from (III) or, less readily, from (II): 


VM -—1 
¢(4, Ma) M-* S$ ¢(4,a+r/M) (M an integer). (VI) 


r= 
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This is particularly valuable as a check, and for computing function values 
for large arguments. It is most conveniently applied when M = 2, 4, 5, 
and 10. 


Computation of the table 

Although in virtue of the difference equation (I) it is sufficient in prin- 
ciple to tabulate the function over a range of the variable comprised 
between any two successive integers, say, a more extended tabulation has 
much greater potential value in practice, since it allows values of the 
function for large arguments to be calculated very expeditiously. A large 
number of decimal places is justified because the functions C{(2n+-1)/2,a' 
can be reached by numerical differentiation: 

dl(s,a)/da sC(s+l1,a). 
In the first place short tables of values were computed for 
a 0-7 (0-05) 0-9 (0-01) 1-1(0-05) 1-3) and a 9-0 (0-1) 10-0, 

using the series (V) and (LV) respectively. The individual terms were 
carried to thirteen places and the sums recorded to twelve. The values | 
of the zeta functions in (V) were taken from Gram’s (4) table, except for 
¢(4) and ¢(3) which had to be known to a higher degree of accuracy. 
C(3) 


calculation gave 


was taken from Gram (5), who gave a 3l-place value, and a new 


C(3) 2-612 375 348 686... . 





From this short series a skeleton table. 
a 0-00 (0-05) 0-90 (0-01) 1-10 (0-05) 3-00 (0-1) 10-0, 


was obtained by applying either the difference equation (I) or the multi- 





plication formula (VI) with M 2 in such a way that every new value 
was arrived at by two independent routes from two or more primitive 
values. The entries were recorded to 11} places and no discrepancies 
greater than 0-3 in the eleventh place were detected. This table was then 
differenced and interpolated to a = 1-00 (0-01) 2-00 (0-02) 5-00 (0-05) 10-00, 
the last half-place being dropped. The eleven-place table showed no 
fluctuations in its differences corresponding to errors as large as three 
units of the last place; however, about half of the entries were checked 
for self-consistency by applying the multiplication formula to groups of 
5, 6, or 11 entries (M 4, 5, or 10 in (VI)). To ten places the results 
were exact. Finally, modified second (central) differences were computed, 
and the whole rounded to ten places. 


It is hoped, therefore, that there are none but occasional rounding-off 
errors in the final table. The fourth central differences are nowhere greater 
than 800 units of the last place, so that the modified second differences 
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alues are sufficient for interpolation throughout the table. It should be noted 
4.5 that the function values are negative and the second differences positive 


in the range covered. 


Since this table was prepared my attention has been drawn to the 


prin existence of an unpublished table by Hensman (6), giving {(s,a) for 
rised s 10 (0-1) 0, a 0(0-1) 2, and (s—1)Z(s, a) for s = 0(0-1)1, a = 0(0-1)2 
1 has to five decimals. Hensman made use of a modification of equation (V) 
the having much improved convergence. 

large — 

ape The table is published by permission of the Chief Scientist, Ministry 


yf Supply : 
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DIRECT CALCULATION OF SMOOTH GUNNERY 
RANGE TABLES 


By J. G. L. MICHEL (National Physical Laboratory, Teddington) 


[Received 19 June 1951] 


SUMMARY 


A technique is developed for inverting non-linear functions of two variables by 
using the linear relations between their derivatives. 


1. Introduction 
In this note a technique is developed for inverting non-linear functions of 
two variables by taking advantage of the linearity of the relationships 
between their differentials. It is illustrated by a specific example from 
external ballistics, but is quite general in its principles. 

For a given projectile and muzzle velocity, the integration of the 
equations of motion gives: 


(1.1) 
where 

x = ground range, 

y = height, 

t = time of flight, 

x = quadrant elevation = dy/dx at t = 0. 

In (1.1) « and y are regarded as dependent variables of the two indepen- 
dent variables ¢, « (a enters into the trajectory calculations merely as a 
parameter, an initial condition). . 

The functions required in gunnery are usually an inversion of those 
obtained by direct integration of the equations of motion, viz. 

t (x,y) | 


, (1.2 
¥ a(a, y) J 


where x and y are regarded as independent variables (arguments). 


2. Rationale of method 


The method consists in obtaining in numerical form values of: 


ot ot 

ae F) ~~. 

OX y const cy x const 
(; ) (: q 

a ’ re 

Oa y const OY} x const 


[Quart. Journ. Mech. and Applied Math., Vol. V Pt. 1 (1952)] 
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simultaneously with the solution of the differential equations, in which x and y 
are again regarded as independent variables. For example, we obtain (as 


a table or graph) (é¢ 0x), 9: 
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Fic. 2. 
From this, by simple quadrature, 
t(x),, 0 Ga dx (2.2) 
1) J \OX}y=0 


is obtained, so giving the time of flight to all points at ground level (see 
Fig. 1). From (€t/0y),, concts © = 0.2 Rei t= N,..., etc., one obtains by 
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quadrature values of time of flight along all vertical lines in Fig. 1, initial 
conditions having already been obtained from (2.2). 

As a check values of time of flight along all horizontal lines in Fig. 2 are 
obtained from (ét/€X),, const by quadrature. 

The values at points of intersection of the lines at Figs. 1 and 2 should, 
of course, be the same. 


Values of a(x, y) are obtained in precisely the same manner by quadra- 


(; “| (- ‘ 
Ov ial cy 


3. Basis of method of calculation 


ture of 


x const 


The differential equations for the trajectory, in which y(¢,«) and éy/ét 
instead of dy/dt, etc., have purposely been used to emphasize the argu- 
ments, are of the form 


9 





072(t, x) [ea(t,a) ey(t, x) ’ 
—_ , y(t, 
ot? J ot ot K\ | | 
a (3.1) 
éty(t, a) éylta) @x(t,a) 4 | | 
ot= ot 7 





with the initial condition (values at t = 0) 


% = yy = 0, 02,/0t = Vcosa, 


Cy, /et VY, sin «, Vy, = muzzle velocity. 

By partially differentiating each of equations (3.1) with respect to « one 
obtains a special case of the equations for the first-order variations from 
a trajectory. It is possible to integrate the resulting equations along a 
trajectory simultaneously with (3.1) to obtain values of éa/éa,- éy/éa along 
a trajectory. 


From these and the solutions of (3.1) it is possible to tabulate values of 


6x/Ot, Cx/Cx, Cy/Ct, Cy/Ex at discrete points where the chosen trajectories cut 
the grid lines of Figs. 1 and 2 (see Fig. 3: cf. also Hartree, Michel, and 
Nicholson (1)). 

From these four numerical values of the partial derivatives when ¢ and « 


are regarded as independent variables, the Implicit Function Theorem 


(Stewart, Advanced Calculus, section 6.32) gives the numerical values of 


the inverted partial derivatives as 


wed ") 
Cox 


(“e2)) _J- — ") | 
cy r const On ‘const 


\ e Vs \ 


J oe x) 


y const Ca ‘const 











DI 


whe 


the 
side 
sim 


par 


col 
pre 


or¢ 








tial 


are 








DIRECT CALCULATION OF SMOOTH GUNNERY RANGE TABLES 1: 


bo 
=I 


(2) (4) 
OX, Y) \at x const ot x const 
O(t, x) re ( 

CX) ¢ const c 2. 


the jacobian, and the numerical values of all quantities on the right-hand 


where J 


side are known. Values of (@«/@2),, const» (€x/@Y)z const are Obtained from 
similar formulae. (A simple derivation of the formulae (3.2) for these 


particular functions is given in the appendix.) 
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4. Conclusion 

Experience suggests that (using a differential analyser of sufficient size) 
values of the partial derivatives (3.2) can be obtained accurately and self 
consistently to the order of three significant figures. Their quadrature can 
produce results consistent (though not, of course, accurate) to a still higher 
order (integration is a smoothing process). 

The comparison of results obtained by the processes exemplified by 
Figs. 1 and 2 at each strip of the grid gives an indication of the over-all 
accuracy. 

The method has been exemplified with reference to the inversion of a 
particular pair of functions of two variables x(t, x) and y(t, x). The principle 
is, of course, perfectly general, and can be applied mutatis mutandis to 
other arguments and their functions. 

Moreover, results are obtained in the form required for gunnery tables 
without inverse interpolation in a double entry table. The principle may 
therefore be of value for digital machines or hand numerical methods, 
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particularly where an integration process such as the Runge-Kutta method, 
in which step size can be varied to give values of the variables on the grid, 
is used. 


APPENDIX 
From (1.1) one obtains 
: 4 on 
ot Cau 


ay) CY « 
oF Se — Ow 


Ca 


oY 
[f y is held constant, dy vanishes 
ba 
0 = 
ot 
The elimination then of da from (A 2) 


(Ox Cy Ox doy\. 
gives - i) 


Ot ca Oa Ct! 


oy | 
| 7 


ot ot 
or lim = { A 
Ox Ox 


vy const 
as in (3.2). 

The expressions for the other partial derivatives are readily obtained by a precisely 
similar process. 
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